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I INTRODUCTION 


This  report  describes  the  MUSAC  II  computer  simulation  model;  the 
model  can  simulate  an  acoustic  classification  decision  process  using 
passive  sonar  information  available  in  scenarios  involving  multiple 
targets.  The  study  objective  was  to  demonstrate  the  MUSAC  II  computer 
program  by  analyzing  a scenario  in  which  a submarine  uses  Lofar  and  Demon 
information  to  classify  a pair  of  ships.  The  purpose  of  the  demons trot ion 
was  to  show  the  capabilities  of  the  computer  program.  This  study,  the 
fourth  of  a series,  completes  the  acoustic  classification  modeling  effort 
that  originated  in  SRI's  evaluation  of  acoustic  countermeasures  concepts 
and  techniques  for  the  Office  of  Naval  Research  (ONR)  in  the  early 
1970s. 

A.  MUSAC  Background 

In  1971,  ONR  initiated  a long-term  research  effort  with  SRI  under 
Contract  N00014-71 -C-0119  to  address  modeling  of  the  acoustic  classifica- 
tion process.  The  general  objective  was  to  explore  alternative  analytic 
methodologies  and  to  recommend  a methodology  to  represent  passive  sonar 
classification  during  a submarine  operation  against  a surface  ship  group. 
More  specifically,  the  methodology  was  to  be  suitable  for  use  in  evaluating 
tactical  deception  concepts,  including  deployment  and  employment  of  acous- 
tic countermeasures  in  the  protection  of  naval  surface  forces.  After  an 
extensive  investigation  into  ways  to  represent  classification,  a methodology 
was  created.  It  was  called  "MUSAC,"  an  acronym  for  "Multiple  Source 
Acoustic  Classification". 

The  objectives  of  a second  ONR  research  task,  started  in  1973  under 
Contract  N00014-71-C-0419,  were  to  recommend  modifications  or  extensions 
for  incorporating  MUSAC  routines  into  existing  large-scale  antisubmarine 
and  antisurface  warfare  engagement  models,  and  to  apply  MUSAC  independently 
for  evaluating  classification  problems.  Although  several  large-scale 


models,  such  as  APSUB,  APSURF,  and  SIM  II,  were  investigated,  the  study 
effort  did  not  succeed  in  incorporating  MUSAC  into  those  models. 

The  basic  product  of  the  second  MUSAC  task  was  a computer  model  that 
was  used  on  two  occasions.  The  first  application  was  the  acoustic 
deception  examples  prepared  for  the  project  report.  The  second  application 
was  in  a Harpoon  targeting  study  and  was  the  first  true  test  of  the  MUSAC 
methodology;  unfortunately,  methodological  difficulties  arose  when  the 
model  was  used  for  convergence  zone  targets.  In  concluding  the  second 
research  effort  on  MUSAC,  a draft  report  was  submitted  to  ONR  for  review. 

At  the  direction  of  the  project  scientific  officer,  the  draft  report  was 
also  reviewed  by  several  Navy  laboratories  and  private  companies.  In 
addition  to  giving  valuable  criticism,  several  respondents  expressed  a 
need  for  a model  like  MUSAC  to  apply  to  acoustic  classification  problems. 

A third  MUSAC  effort  was  started  in  1975  under  Contract  N0014-76-C- 
0166.  The  objective  was  to  enhance  the  MUSAC  methodology  by  incorporating 
the  suggestions  received  from  its  review,  and  to  solve  several  problems 
encountered  during  the  application  of  the  model.  The  project  objectives 
were  to  generalize  and  modify  the  existing  computer  programs,  to  analyze 
and  resolve  methodologic  questions,  and  to  revise  and  expand  MUSAC 
documentation.  The  third  MUSAC  effort  was  not  completed,  because  the 
study  direction  was  changed  by  ONR  about  halfway  through  the  project.  As 
part  of  an  ONR  reorganization,  the  MUSAC  project  was  transferred  from 
Code  431  to  Code  230.  After  the  new  scientific  officer  evaluated  the 
methodology's  potential  for  his  Fleet-oriented  program,  the  project  was 
redirected  to  a tactical  development  task  totally  unrelated  to  the  goals 
of  MUSAC.  Even  with  only  half  the  initial  funds,  good  progress  was  made 
toward  achieving  the  objectives. 

The  MUSAC  part  of  the  third  project  produced  a completely  revised 
and  documented  methodology.  The  MUSAC  methodology  was  restructured 
extensively  enough  that  it  was  called  "MUSAC  II"  to  distinguish  it  from 
the  early  methodology.  Major  revisions  included  a new  formulation  of  a 
multifeature  sonar  detection  model,  different  likelihood  calculations, 
and  a more  generalized  decision-making  procedure.  Although  the  MUSAC  II 
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methodology  was  well -documented,  there  were  not  enough  funds  to  implement 
the  methodology  by  developing  the  computer  program  and  demonstrating  its 
capabilities.  Revising  the  original  computer  model  was  not  feasible 
because  of  the  many  basic  changes  in  the  methodology.  Thus,  MUSAC  II 
required  a new  computer  program. 

In  1976,  SRI  International  supported  an  IR6cD  project  to  develop  a 
MUSAC  II  computer  program.  The  program  was  coded,  keypunched,  and 
corrected  for  compilation  errors;  but  it  was  not  demonstrated  with  specific 
input  parameters  and  functions. 

Demonstration  of  the  computer  program  was  the  goal  of  the  present 
MUSAC  project;  the  program  and  its  capabilities  are  documented  in  this 
report . 


B.  MUSAC  II  Application 

The  MUSAC  II  model  may  be  (1)  used  as  a component  of  a larger 
Monte  Carlo  acoustic  warfare  engagement  simulation,  or  (2)  used  by  itself 
as  an  analytical  tool.  In  the  first  application,  the  model  can  provide 
classification  decisions  for  the  selection  of  tactics  in  dynamic  engage- 
ments in  which  the  simulated  kinematics  are  not  predetermined.  In  the 
second  application,  the  model  can  predict  classification  probabilities 
and  perform  sensitivity  analyses  of  acoustic  parameters.  MUSAC  II  employs 
the  standard  acoustic  parameters  of  classical  sonar  detection  theory.  By 
using  a physical-based  approach,  the  model  represents  the  inherent  classi- 
fication capability  of  a sonar  system,  particularly  the  sensitivity  of 
the  classification  decisions  to  signal-to-noise  ratios. 

Large,  detailed  computer  simulations  of  the  future  are  expected  to 
be  an  active  research  medium  because  of  constant  improvements  in  computer 
capabilities  and  speed.  Since  the  MUSAC  II  approach  is  much  more 
sophisticated  than  the  ad  hoc  acoustic  classification  algorithms  currently 
used  in  engagement  simulations,  the  model  will  be  useful  in  future  sonar 
systems  analyses  using  computer  simulation. 

1.  J.  R.  Olmstead  and  T.  R.  Elfers,  "MUSAC  II,  A Method  for  Modeling 
Passive  Sonar  Classification  in  a Multiple  Target  Environment," 
NWRC-TN-62,  SRI  International,  Menlo  Park,  California  (February  1976), 
UNCLASSIFIED,  AD-A028-936. 

3 


C.  MUSAC  II  Methodology  and  Computer  Program 


The  MUSAC  II  methodology  is  a mathematical  representation  of  passive 
sonar  classification,  and  the  principal  attribute  of  the  methodology  is 
its  multiple-target  capability.  Almost  without  exception,  other  models 
allow  for  only  one  target  at  a time.  The  methodology  is  based  on  the 
detection  of  acoustic  features;  in  this  way,  spectral  and  spatial 
acoustic  information  is  modeled  so  that  the  sonar  system's  bearing  and 
frequency  resolution  influences  the  classification  outcome.  The  acoustic 
features  are  defined  by  the  analyst;  the  features  can  be  narrowband, 
broadband,  or  modulated  broadband  classification  clues  (for  example, 

Lofar  or  Demon  lines).  The  acoustic  features  are  represented  by 
Bernoulli  random  variables  so  that  the  stochastic  structure  of  the  model 
provides  for  realistic  random  variations  of  acoustic  data.  A dynamic 
encounter  is  represented  by  a time-step  simulation.  The  MUSAC  II 
methodology  is  structured  for  sequential  decision-making  by  the  update 
of  classification  information  and  the  change  in  kinematic  variables  over 
time.  From  Monte  Carlo  replications  of  the  encounter,  the  probability 
of  making  selected  tactical  and  classification  decisions  can  be  estimated. 

The  MUSAC  II  methodology  uses  a Bayesian  decision-making  approach. 

The  analyst  first  formulates  a set  of  multiple-target  hypothesis  that 
will  be  used  in  the  engagement  simulation.  The  probability  of  detecting 
specified  acoustic  features  is  calculated  at  each  time  step,  for  each 
sonar  look  angle,  and  for  each  hypothesis  (the  true  target  configuration 
is  one  of  the  hypotheses).  These  detection  probabilities  are  then  used, 
in  conjunction  with  the  observed  random  features,  to  calculate  the 
likelihood  that  the  data  would  be  observed  if  the  hypothesis  were  true. 

The  likelihoods  and  the  prior  probabilities  are  then  combined  to  produce 
the  posterior  probability  that  the  Ith  hypothesis  is  true,  given  the 
observed  data.  The  analyst  defines  tactical  or  classification  decisions 
that  are  to  be  simulated,  the  value  of  making  the  decision  when  each 
hypothesis  is  true,  and  value  thresholds.  With  this  decision  structure, 
MUSAC  II  determines  if  a decision  is  to  be  made  at  the  present  time 
step;  if  not,  another  time  step  is  simulated  and  more  data  collected. 


The  computer  p og’am  is  coded  in  Fortran  Extended  ^or  SRI's  CDC 
6400  computer.  The  source  deck  is  approximately  1,000  cards;  and 
with  the  currently  programmed  array  dimensions,  the  program  requires 
about  32,800  words  of  memory.  The  running  time  is  directly  proportional 
to  the  number  of  runs,  number  of  replications,  and  number  of  time  steps; 
the  running  time  is  influenced  to  a lesser  extent  by  the  number  of 
target  tracks,  number  of  hypotheses,  and  number  of  features.  For  example, 
with  the  parameters: 

3 runs 

20  replications 
10  time  steps 
2 target  tracks 
9 hypotheses 
7 features, 


the  running  time  was  121  seconds. 


II  SINGLE -TARGET  EXAMPLE 


This  chapter  describes  how  the  MUSAC  II  model  can  simulate  a single- 
target scenario;  the  next  chapter  deals  with  the  double  target  case. 

The  input  parameters  are  first  explained,  then  the  model  output  is 
discussed,  and  fina1 ly  variations  on  the  single-target  scenario  are 
investigated  to  show  the  sensitivity  of  classification  probabilities  to 
various  input  parameters.  Appendix  A contains  a listing  of  the  computer 
model.  Appendix  B contains  (1)  a listing  of  Subroutine  INPUT,  which 
shows  the  input  parameter  values  used  in  the  example;  and  (2)  a listing 
of  the  output,  which  shows  the  results  of  the  single-target  simulation. 

The  single-target  scenario  involves  classification  of  a single 
surface  ship  by  a submarine  using  Demon  information  from  a hull -mounted 
array  and  Lofar  information  from  a towed  array.  The  initial  range 
separation  between  the  target  ship  and  the  submarine  is  about  40  nmi,  and 
the  two  units  approach  each  other  on  a near  collision  course.  The  engage- 
ment lasts  for  about  1.5  hours  and  ends  when  the  units  are  6 nmi  apart. 
During  the  1.5-hour  period,  the  submarine  must  classify  the  target  as  one 
of  three  classes  of  surface  ships. 

A.  Input  Parameters  and  Functions 

The  model  requires  50  parameters  and  six  user-defined  functions. 

Many  of  the  parameters  are  multivariate,  meaning  that  they  are  arrays 
using  one  or  two  subscripts.  Table  1 lists  the  input  parameters  and 
functions.  The  input  list  is  subdivided  into  seven  categories.  The 
following  sections  discuss  what  the  parameters  mean,  how  they  are  used, 
and  what  value  they  assume  in  the  single-target  example  computation. 

The  input  parameters  for  the  model  are  contained  in  Subroutine 
INPUT  in  the  form  of  DATA  statements;  this  method  of  inputting  parameters 
was  chosen  because  of  the  versatility  of  Fortran  DATA  statements,  even 
though  a small  price  is  paid  in  recompiling  INPUT  each  time  the  program 
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Table  1 

INPUT  PARAMETERS  AND  FUNCTIONS 

1.  Scenario  Parameters 


GV 

Random  number  generative  value 

NREP 

Number  of  replications 

NMAX 

Number  of  time  steps 

TS 

Time  step  duration  (min) 

IMAX 

Number  of  hypotheses 

MMAX 

Number  of  tracks 

KT(I,M) 

Hypothesized  target  type 

IR 

True  hypothesis  number 

PRIOR(I) 

Prior  probability  of  hypothesis 

KDMAX 

Number  of  tactical  decision  alternatives 

VAL(I,KD) 

Value  of  decision 

FVAL 

Decision  threshold  of  test  ratio 

XT(M) 

YT(M) 

CT(N,M) 

ST(N,M) 

DT(N,M) 


Target  initial  x-position  (nmi) 
Target  initial  y-position  (nmi) 
Target  course  (deg) 

Target  speed  (kt) 

Target  depth  (ft) 


Target  Classification  Parameters 


JMAX  Number  of  features 

NF(J)  Feature  off/on  (0  = off  1 = on) 

KF(J)  Feature  type  (1  = Lofar  2 = BBand  3 = Demon) 

FRQ(J)  Center  frequency  (Hz) 

PLL(J,KT)  Lofar  line  level  (dB  //  uPa  at  1 yd)  also 

PLL(J,KT)  Demon  modulation  level  (dB) 

PBB(NB.KT)  Broadband  source  spectrum  (dB  //  uPa  /Hz  at  1 yd) 

FQB(NB,KT)  Frequency  points  for  PBB  (Hz) 


Table  1 (Continued) 


4.  Sonar  Array  Parameters 


LAMAX 

Number  of  arrays 

LA(J) 

Array  number 

KA(LA) 

Array  type  (1  = circle,  2 = line) 

DH(LA) 

Horizontal  array  size  (ft) 

DV(LA) 

Vertical  array  size  (ft) 

HN(LA) 

Horizontal  number  of  hydrophones 

VN(LA) 

Vertical  number  of  hydrophones 

SL(LA) 

Sidelobe  level  (ratio) 

XO(LA) 

Array  initial  x-position  (nmi) 

YO(LA) 

Array  initial  y-position  (nmi) 

C0(N, LA) 

Array  course  (deg) 

S0(N,LA) 

Array  speed  (kt) 

D0(N, LA) 

Array  depth  (ft) 

NA(N,LA) 

Array  off/on  (0  = off,  1 = on) 

2 

PNN(NF,LA) 

Broadband  noise  outside  array  (dB 

//  u Pa  /Hz) 

FQN(NF.LA) 

Frequency  points  for  PNN  (Hz) 

5.  Signal  Processing  Parameters 


LP(J) 

Processor 

number 

WTH(LP) 

Bandwidth 

(Hz) 

SCR(LP) 

Scan  rate 

(nurr^er/min) 

FCS(LP) 

Number  of 

feature  cells  per  scan 

TOT(LP) 

Total  time 

; of  feature  integration  (min) 

DET(LP) 

Detection 

threshold  (number  of  sigmas) 

6.  Acoustic  Fluctuation  Parameters 


MIX(KR)  Mixing  constant  (0  = gauss,  1 = jump) 

SDV(KR)  Standard  deviation  of  random  process  (dB) 

TAU(KR)  Relaxation  time  of  random  process  (min) 

KR:  1 = Lofar  4 = Noise 

2 = BBand  5 = PLoss 

3 = Demon 


7.  User-Defined  Functions 


FLL(J , KT, ST, DT, ASP) 
FBB(KT,FRQ,ST,DT,ASP) 
FNN ( LA , FRQ , S 0 , DO , ANG ) 
FAZ(FRQ, RNG.DO, DT) 
FBM(J,FRQ, BRG.ANG) 
FSL(LA.FRQ) 


Lofar  and  Demon  source  level 
Broadband  source  spectrum  levei. 

Noise  spectrum  level  at  array  output 
Propagation  loss 
Beam  pattern  ratio 
Reference  beam  pattern  ratio 
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is  run.  The  input  functions  are  coded  as  individual  Fortran  functions, 
thus  allowing  for  easy  changes. 

1 . Scenario  Parameters 

GV  is  the  generative  value  used  in  the  library  subroutine 
RANSET(GV)  to  start  a sequence  of  random  numbers  through  use  of  the 
function  RANF.  The  value  of  GV  can  be  set  to  any  large  positive  number, 
for  example  GV  = 583.  The  model  is  programmed  so  that  the  random  number 
sequence  in  one  replication  is  independent  of  the  sequence  in  another 
replication.  However,  the  random  number  sequence  is  not  independent  from 
run  to  run.  By  using  the  parameter  GV,  the  same  sequence  of  random 
numbers  is  used  from  one  run  to  another,  providing  that  each  run  calls 
the  random  number  generator  the  same  number  of  times.  The  purpose  of 
repeating  the  sequence  is  to  allow  for  parameter  sensitivity  analyses 
that  reflect  only  the  variation  of  the  parameters,  not  the  randomness 
of  the  model. 

NREP  is  the  number  of  Monte  Carlo  replications  used  to  compute 
statistics  of  the  engagement.  The  single-target  example  computation  uses 
20  replications,  although  10  times  as  many  would  be  preferred.  Since  the 
purpose  of  the  study  was  model  demonstration  rather  than  model  accuracy, 
a small  number  of  replications  was  adequate;  a small  value  of  NREP  allowed 
for  more  model  demonstration  runs  because  the  cost  per  run  was  less. 

NMAX  is  the  number  of  time  steps  in  the  engagement,  and  TS  is 
the  duration  of  each  time  step.  The  example  computation  used  10  time 
steps  of  10  minutes  each,  for  a total  engagement  time  of  100  minutes. 
Fairly  large  time  steps  were  used  to  reduce  computer  costs.  The  duration 
should  be  set  so  that  (1)  the  target  will  not  "jump  over"  phenomena  such 
as  convergence  zones,  (2)  the  total  integration  time  is  not  shorter  than 
the  time  step,  and  (3)  the  geometry  will  not  change  significantly  between 
time  steps.  The  last  point  relates  to  the  problem  of  using  the  geometry 
at  a point  in  time  as  an  approximation  of  an  average  geometry  over  the 
time  step.  The  model  assumes  that  the  geometry  at  the  end  of  the  time 
step  is  adequate  for  simulating  the  integrative  processes  over  the  whole 
time  step. 
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IMAX  is  the  total  number  of  hypothetical  target  configurations. 
For  the  single-target  example  there  are  three  hypotheses:  (1)  the  target 
is  Type  1,  (2)  the  target  is  Type  2,  and  (3)  the  target  is  Type  3.  For 
multiple  target  cases,  the  track  number  must  also  be  specified  in  the 
definition  of  the  hypothesis.  For  example,  a hypothesis  might  read: 

The  target  on  Track  2 is  Type  3. 

MMAX  is  the  total  number  of  target  tracks  in  the  simulation; 

MMAX  = 1 for  the  single-target  example.  Tracks  are  thought  of  as 
entities  unto  themselves;  when  different  types  of  targets  are  placed  on 
the  tracks,  different  hypotheses  are  generated.  The  computer  model  is 
simplified  by  allowing  only  the  tracks  of  the  true  target  configuration 
to  be  used  in  forming  the  hypotheses.  Thus  there  is  an  underlying 
assumption  that  the  geometry  of  the  engagement  situation  is  known.  This 
simulated  knowledge  may  have  an  actual  basis  as  objective  knowledge  (a 
tracking  solution)  or  subjective  knowledge  (long-range  targets  imply 
low  signals  and  frequency  attenuation).  The  simulation  methodology  must 
use  reasonable  geometries  for  computations;  since  the  computations  were 
already  overburdened  with  replication,  time-step,  and  feature  calculations, 
the  inclusion  of  pseudo-tracks  was  not  attempted. 

KT( I,M)  is  an  array  that  defines,  for  the  I-th  hypothesis,  what 
type  of  target  is  on  the  M-th  track.  In  the  single-target  example 
KT(I,1)  = 1,2,3  for  I = 1,2,3.  Thus  for  the  first  hypothesis,  target 
Type  1 is  on  Track  1;  for  the  second  hypothesis,  target  Type  2 is  on 
Track  1;  and  so  on.  As  currently  programmed,  up  to  10  hypotheses  can  be 
defined  over  two  tracks;  however,  these  array  dimensions  can  easily  be 
changed.  There  is  no  computer  restriction  on  the  number  of  target  types, 
since  "type"  is  the  value  of  the  array,  but  many  types  implies  many 
hypotheses,  so  in  effect  the  number  of  target  types  is  limited  to  the 
number  of  hypotheses  (or  less,  in  the  case  of  multiple  target  configura- 
tions ) . 

IR  is  the  hypothesis  that  represents  the  real  configuration.  In 
the  single-target  example  IR  = 1 for  the  first  run;  thus  a target  of 
Type  1 is  actually  present,  and  it  may  be  classified  as  Type  1,  2,  or  3 
by  choosing  Hypothesis  1,  2,  or  3.  The  single-target  example  makes 
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three  separate  runs  for  IR  = 1,2,3  so  that  a 3-by-3  classification  matrix 
can  be  formed.  The  MUSAC  II  methodology  does  not  require  that  the  real 
configuration  be  present  as  one  of  the  hypothetical  configurations; 
however,  the  computer  model  was  simplified  by  designating  one  of  the 
hypotheses  as  true.  Also,  the  interpretation  of  "correct"  classification 
is  more  clear  when  one  of  the  hypotheses  is  true. 

PRIOR( I)  is  proportional  to  the  a priori  probability  that  the 
I-th  hypothesis  is  true.  In  the  single-target  example,  the  priors  are 
equal  to  0.1  for  all  three  hypotheses  (the  priors  do  not  have  to  add  to 
1.0,  since  they  are  used  as  weights  in  the  calculations).  A priori 
knowledge,  such  as  order-of-battle  estimates  or  historical  track  records, 
can  be  modeled  by  appropriately  chosen  priors. 

KDMAX  is  the  total  number  of  tactical  decisions  alternatives. 

In  the  first  part  of  the  example  calculation,  KDMAX  = 3 to  correspond 
with  the  three  possible  target  types.  For  this  case  the  "tactical 
alternatives"  are  decisions  to  classify  the  target  as  Type  1,  2,  or  3. 
Later  examples  have  KDMAX  = 2,  and  the  tactical  alternatives  are  (1)  to 
attack,  or  (2)  to  break  off  the  approach.  As  currently  programmed,  up 
to  10  decision  alternatives  can  be  defined. 

VAL( I,KD)  is  an  array  that  specifies  the  value  of  decision 
alternative  KD,  given  that  the  I-th  hypothesis  is  true.  When  the 
decision  alternatives  are  to  classify  the  target,  the  VAL(I,KD)  matrix  is 


KD 

1 

0 

0 

I 

0 

1 

0 

0 

0 

1 

In  other  words,  a value  of  1 is  assigned  to  a correct  classification  and 
a value  of  zero  to  an  incorrect  classification  decision.  When  the 
decision  is  to  attack  or  not,  the  VAL(I,KD)  matrix  is: 
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KD 

300 

-500 

I 

-100 

100 

-500 

100 

where  the  first  column  is  the  "attack"  decision.  The  values  imply  that  a 
high  penalty  is  paid  for  attacking  Type  2 and  3 targets;  however,  a high 
penalty  is  also  paid  for  not  attacking  a Type  1 target. 

FVAL  is  the  threshold  against  which  a computed  ratio  is  tested. 
The  ratio  is  the  "maximum  expected  value"  divided  by  the  "expected 
maximum  value."  "Value"  refers  to  the  decision  alternative  values 
VAL(I,KD),  and  "expected"  means  that  the  values  are  averaged  by  using 
the  posterior  probability  distribution  POST(I).  The  test  ratio  is 
defined  as: 

Max  [Z  POST(I)  VAL(I,KD)] 

F i , 

Z POST(I)  MaxtVAL(I.KD)] 

where  the  maximizing  operation  is  over  the  decision  alternatives  KD.  The 
F ratio  is  between  0 and  1,  and  for  the  example  calculation  FVAL  = 0.95. 
When  F is  less  than  FVAL,  the  decision  is  deferred  and  more  information 
is  collected  by  letting  the  model  advance  another  time  step;  when  F is 
equal  to  or  greater  than  FVAL,  a decision  is  made.  The  chosen  decision 
alternative  is  the  one  that  corresponds  to  the  maximum  expected  value: 

Select  alternative  KD*  such  that: 

■ 

E(KD* ) = Max  E(KD)  . 

E(KD)  is  defined  as  the  expected  value: 

E(KD)  = I POST(I)  VAL(I,KD)  . 

The  above  decision-making  method  is  slightly  different  from  that  described 
in  the  MUSAC  II  methodology  report;^  there,  the  difference  of  numerator 
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and  demonimator  was  tested  instead  of  the  ratio.  One  way  works  as  well 
as  the  other,  and  the  present  method  has  the  advantage  of  using  a 
dimensionless  input  parameter,  FVAL,  that  does  not  have  to  be  changed 
when  the  VAL  matrix  is  changed. 

2.  Target  Track  Parameters 

XT(M)  and  YT(M)  specify  the  initial  position  of  target  Track  M. 
The  trajectory  through  time  is  defined  by  CT(N.M)  and  ST(N,M) , the  course 
and  speed  values  during  the  N-th  time  step.  The  model  accepts  target 
depth,  DT(N.M).  as  an  input  parameter,  but  does  not  currently  use  it  in 
any  calculations.  The  primary  use  of  depth  would  be  in  the  calculation 
of  propagation  loss;  however,  a depth-independent  propagation  model  is 
used  for  model  demonstration.  The  single-target  example  uses  a straight- 
running target  and  has  the  following  input  parameters: 

XT(1 ) = 0 nmi 

YT(1)  = 40  nmi 

CT(l,N)  = 155  deg  N = 1,10 

ST(1,N)  = 25  kt  N = 1,10 

3 . Target  Classification  Parameters 

JMAX  is  the  total  number  of  classification  features  that  are 
used  to  classify  targets.  The  example  uses  seven  features:  three  Lofar 
lines,  and  two  Demon  lines  in  two  bands  (a  total  of  four  Demon  features). 
As  currently  programmed,  up  to  12  features  may  be  used.  Any  combination 
of  Lofar  lines,  Demon  lines  in  various  bands,  and  broadband  noise  may 
be  defined. 

NF(J ) is  an  array  that  singles  out  classification  features  to 
be  used  or  ignored.  When  NF(J)  = 1 the  J-th  feature  is  used,  and  when 
NF(J)  = 0 the  J-th  feature  is  treated  as  though  it  did  not  exist.  NF 
used  to  study  the  importance  of  individual  features  by  turning  them  off 
and  then  on  in  successive  runs. 
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J is : 


KF( J ) tells  the  model  what  kind  of  classification  feature 


KF  = 1 Lofar  line 

2 Broadband  noise 

3 Demon  line 


In  the  example  run: 


KF(J)  = 1 for  J = 1,2,3 
KF( J)  = 3 for  J = 4, 5, 6, 7 , 

Broadband  noise  is  included  as  a feature  in  case  the  analyst  desires 
the  "loudness"  of  the  target  to  convey  classification  information. 

FRQ( J ) is  the  center  frequency  of  the  J-th  classification 
feature.  For  example,  the  three  Lofar  lines  are  at  50  Hz,  100  Hz,  and 
400  Hz,  and  the  two  Demon  lines  are  in  two  bands  centered  at  2,828  Hz 
and  5,656  Hz  (the  geometric  mean  of  the  2 to  4 kHz  band  and  the  4 to  8 kHz 
band).  The  frequency  parameter  is  used  primarily  in  the  attenuation 
calculation  associated  with  propagation  loss  and  in  the  beam  pattern 
calculation. 

PLL( J,KT)  is  the  line  level  of  classification  feature,  J,  for 
target  type,  KT.  For  example,  the  Lofar  line  levels  for  Type  1,  2,  and 
3 targets  are: 


Feature 

J 

Lofar 

Line 

Target  Line  Levels  (dB) 
KT  = 1 2 3 

1 

50  Hz 

155 

155 

0 

2 

100  Hz 

0 

150 

150 

3 

400  Hz 

150 

0 

150 

2 

where  the  Lofar  line  levels  are  in  units  of  dB  relative  to  1 pPa  at 
1 yard.  The  very  small  value  of  0 dB  represents  a missing  line;  for 
example,  Type  1 target  has  lines  at  50  Hz  and  400  Hz,  but  none  at 
100  Hz.  The  Demon  line  levels  are  also  defined  in  the  PLL  matrix  for 
Type  1,  2,  and  3 targets;  for  example: 


Feature 

J 

Demon 

Line 

Band 

(kHz) 

Target  Modulation  Levels  (dB) 

KT  = 1 2 3 

4 

A 

2 to  4 

- 3 

0 

CN 

1 

-20 

5 

B 

2 to  4 

- 3 

- 3 

0 

CN 

1 

6 

A 

4 to  8 

- 3 

-20 

- 2 

7 

B 

4 to  8 

-20 

- 3 

- 2 

where  the  Demon 

line  levels  are 

decibel  values 

of  the 

square  of  the 

modulation  index 

. The  large  negative  modulation  levels  represent 

missing  Demon  lines. 

PBB(NB.KT)  is  the  broadband  source  spectrum  level  at  frequency 
points,  FOB ( NB . KT ) for  target  type,  KT.  The  spectrum  is  described  by  a 
piecewise  linear  function  with  breakpoints  NB;  a maximum  of  six  points 
may  be  specified.  In  the  example  calculation,  the  spectrum  is  assumed 
to  be  identical  for  all  three  target  types  (for  KT  = 1,2,3)  and  is 
described  by: 

NB  FQB(NB,KT)  PBB(NB,KT) 


1 

10 

Hz 

153 

dB 

2 

100 

Hz 

155 

dB 

3 

1,000 

Hz 

145 

dB 

4 

10,000 

Hz 

125 

dB 

2 

where  the  spectrum  is  in  units  of  dB  relative  to  1 (jPa  /Hz  at  1 yard. 
Figure  1 shows  the  broadband  spectrum  levels  and  Lofar  lines  for  the 
Type  1 target. 

4.  Sonar  Array  Parameters 

LAMAX  is  the  number  of  sonar  arrays  in  the  model . As  currently 
programmed,  there  can  be  a maximum  of  two  arrays,  and  the  example  calcula- 
tion uses  both  of  them  (LAMAX  = 2). 
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FIGURE  1 TYPE  1 TARGET  SPECTRUM 

LA( J ) specifies  from  which  array  Feature  J is  derived.  In  the 

example: 

LA(J)  = 2 for  J = 1,2,3 
LA(J)  = 1 for  J = 4, 5, 6, 7 . 

Thus  the  Lofar  features  come  from  Array  2,  and  the  Demon  features  come 
from  Array  1 . 

KA(LA)  specifies  the  type  of  Array  LA.  There  are  two  types 
currently  programmed:  KA  = 1 (circle)  represents  a cylindrical  or 
spherical  array,  and  KA  = 2 (line)  represents  a towed  array.  In  the 
example  the  first  array  is  circular  and  the  second  is  linear: 

KA( 1 ) = 1 
KA( 2 ) = 2 . 

Thus  the  Lofar  features  are  derived  from  the  towed  array,  and  the  Demon 
features  from  the  hull-mounted  cylindrical  array.  The  difference  between 
the  two  types  of  arrays  is  in  the  way  the  beam  pattern  is  calculated,  as 
described  in  the  section  on  Function  FBM. 
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DH(LA)  is  the  horizontal  dimension,  and  DV(LA)  is  the 
vertical  dimension  of  Array  LA.  In  the  example  calculations: 

LA  DH(LA)  DV(LA) 


60  ft 


Thus  the  cylindrical  array  is  7 ft  in  diameter  and  5 ft  high,  and  the 
towed  array  is  60  ft  long  (the  vertical  dimension  for  a line  array  is 
ignored ) . 

HN(LA)  is  the  number  of  hydrophones  in  the  horizontal  direction, 
and  VN(LA)  is  the  number  of  hydrophones  in  the  vertical  direction  for 
Array  LA.  In  the  example: 

LA  HN(LA)  VN(LA) 


Thus,  the  cylindrical  array  has  15  x 11  hydrophones  that  can  receive 
in  any  one  direction,  and  the  towed  array  has  30  hydrophones.  The 
number-of -hydrophones  parameter  is  used  in  function  FNN  to  calculate  the 
directivity  index. 

SL(LA)  is  the  sidelode  level  of  Array  LA;  it  is  input  as  a 
ratio  rather  than  as  a negative  decibel  value.  The  example  uses 
SL  = 0.01  for  both  arrays;  this  corresponds  to  a maximum  sidelode  level 
of  -20  dB.  The  function  that  calculates  the  beam  pattern  uses  SL  as  a 
limiting  value  so  that  the  beam  response  is  at  least  20  dB  below  the 
main  beam  response  whenever  the  signal  arrives  outside  the  main  beam. 

XQ(LA)  and  YO(LA)  specify  the  initial  position  of  Array  LA.  The 
array's  track  is  defined  by  CO(N.LA)  and  S0(N, LA),  the  course  and  speed 
of  Array  LA  during  time  step  N.  The  model  accpets  the  array  depth 
DO(N.LA)  as  an  input  parameter,  but  does  not  currently  use  it  in  any 
calculation.  As  with  target  depth,  array  depth  could  be  used  in  propaga- 
tion-loss calculations.  The  tracks  of  the  two  arrays  in  the  model 
demonstration  are  given  by: 


Parameter 

LA  = 1 

LA 

= 2 

XO(LA) 

0 nmi 

-0.5 

nmi 

YO(LA) 

0 nmi 

0 

nmi 

CO(N, LA) 

90  deg 

90 

deg 

so(n, LA) 

10  kt 

10 

kt 

The  towed  array  follows  1/2  nmi  behind  the  cylindrical  array;  both 
move  east  at  10  kt. 

NA(N,LA)  specifies  whether  or  not  Array  LA  is  operating  during 
time  step  N:  NA  = 0 means  that  the  array  is  off,  and  NA  = 1 means  that 
the  array  is  on.  The  NA(N,LA)  matrix  can  be  used  to  turn  off  a towed 
array  during  a course  change,  or  to  investigate  classification  informa- 
tion on  one  array  versus  that  on  another.  The  example  calculation  leaves 
both  arrays  on  during  the  entire  encounter. 

PNN(NF.LA)  is  the  broadband  noise  spectrum  level  outside 
Array  LA  at  frequency  points  FQN(NF,LA) . The  spectrum  is  described 
by  a piecewise  linear  function  with  breakpoints  NF;  a maximum  of  six 
points  may  be  specified.  In  the  example  calculation,  the  noise  spectra 
outside  the  two  arrays  are  described  by: 


FQN(NF.l)  PNN(NF.l) 


FQN(NF, 2 ) PNN(NF, 2 ) 


10  Hz  75  dB 

1,000  Hz  65  dB 

10,000  Hz  45  dB 

2 

where  the  noise  spectra  are  in  units  of  dB  relative  to  1 ^Pa  /Hz. 

5 . Signal  Processing  Parameters 

LP( J)  is  the  signal  processor  that  produces  Feature  J;  as 
currently  programmed,  LP  may  range  from  1 to  5.  The  example  uses  three 
processors  to  correspond  to  the  different  signal  processing  parameters 
for  Lofar  signals  (LP  = 1),  and  Demon  signals  in  two  bands  (LP  = 2,3). 

WTH(LP)  is  the  bandwidth  of  processor  LP.  The  bandwidth  should 
be  set  to  the  natural  width  of  the  signal  being  processed;  the  model  does 
not  have  a mechanism  to  decrease  the  Lofar  line  level  when  the  bandwidth 
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10  Hz  120  dB 

300  Hz  65  dB 

1,000  Hz  65  dB 

10,000  Hz  45  dB 


is  too  narrow.  A bandwidth  larger  than  necessary  will  allow  more  noise 
into  the  system  and  degrade  detection  performance.  The  example  Lofar 
lines  are  assumed  to  be  1 Hz  wide,  and  the  Demon  bands  are  2 kHz  and  4 kHz 
wide : 


LP  WITH(LP) 

1 1 Hz 

2 2,000  Hz 

3 4,000  Hz  . 

SCR(LP)  is  the  scan  rate  of  the  processor;  it  is  the  number  of 
times  per  minute  that  a particular  frequency  or  bearing  region  is 
processed.  For  example,  the  Lofar  lines  are  assumed  to  lie  in  a 0 to 
600  Hz  frequency  region  that  is  scanned  once  per  minute,  and  the  Demon 
lines  are  assumed  to  lie  in  a 0 to  60  Hz  modulated  frequency  region  that 
is  also  scanned  once  per  minute: 

LP  SCR(LP) 

1 1 scan/min 

2 1 scan/min 

3 1 scan/min 

The  number  of  times  a feature  is  scanned  in  one  time  step  is  calculated 
by: 


NUM  = SCR  * TS  . 


In  the  example,  there  are  10  independent  observations  (scans) 
during  each  time  step. 

FSC(LP)  is  the  number  of  feature  cells  in  one  scan  for  Processor 
LP.  In  the  Lofar  case,  one  scan  covers  600  Hz  and  the  cell  size  is  1 Hz, 
therefore,  FCS  = 600.  For  the  Demon  case,  it  is  assumed  that  the 
resolution  is  1 Hz  of  modulated  frequency,  and  a 60-Hz  scan  implies  60 
feature  cells: 


LP  FSC(LP) 

1 600  cells 

2 60  cells 

3 60  cells 

By  using  the  scan  rate  and  number  of  features  per  scan,  the  averaging 

time  per  scan  for  a given  feature  is: 

TAV  = 60/ (SCR  * FCS ) seconds. 

This  is  the  averaging  time  in  the  time-bandwidth  product  used  in  calcula- 
ting the  standard  deviation  of  the  output  of  the  signal  processor.  In 
the  example,  the  integration  time  per  scan  for  a Lofar  feature  is  0.1  s, 
and  the  integration  time  per  scan  for  Demon  feature  is  1 s. 

TOT(LP)  is  the  total  amount  of  integration  time.  For  example, 
the  Lofar  signals  are  assumed  to  be  traced  on  a moving  paper  recorder  and 
about  30  min  worth  of  visual  integration  is  available.  The  Demon  signals 
are  assumed  to  have  1 min  of  integration: 

LP  TOT(LP) 

1 30  min 

2 10  min 

3 10  min 

The  model  allows  for  information  to  be  accumulated  over  several  time 
steps.  The  number  of  time  steps  that  are  remembered  past  the  current 
time  step  is: 


MEM  = (TOT  - TS )/TS 

The  memory  for  Lofar  features  is  two  time  steps,  and  the  Demon  features 
have  zero  memory. 

The  signal  processor  input  parameters  are  in  terms  of  scan 
rate,  features  per  scan,  and  total  integration  time  because  these  are 
systems  parameters.  The  MUSAC  II  methodology,  however,  uses  averaging 
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time  per  observation,  number  of  independent  observations  per  time  step, 
and  number  of  time  steps  to  be  remembered;  these  are  the  NUM,  TAV,  and  MEM 
parameters  calculated  using  the  model  inputs  SCR,  FCS,  and  TOT. 

DET(LP)  the  number  of  sigmas  above  the  mean  reference  output 
where  the  feature  detection  threshold  is  set.  The  example  calculation 
uses  DET  = 2 sigmas  for  all  three  processors.  The  processor  output  of  the 
data  channel  is  the  simulated  average  square  pressure.  The  random  output 
is  drawn  from  a normal  distribution  and  compared  to  the  threshold  value; 
if  larger  than  the  threshold,  then  the  feature  is  detected.  The  mean 
and  sigma  values  used  to  set  the  threshold  value  are  based  on  the  statis- 
tics of  the  reference  channel;  these  reference  statistics  are  usually 
different  from  the  statistics  of  the  data  channel.  The  feature  detection 
model  is  described  in  the  MUSAC  II  methodology  report.^ 

6 . Acoustic  Fluctuation  Parameters 

MIX(KR)  is  a parameter  that  determines  the  amount  of  Gaussian 
fluctuation  versus  lambda-sigma  jump  fluctuations  for  acoustic 
phenomenon  KR: 

KR  Fluctuation  in 

1 Lofar  line  level 

2 Broadband  source  spectrum  level 

3 Demon  modulation  level 

4 Broadband  noise  spectrum  level 

5 Propagation  loss 

These  five  phenomena  are  simulated  by  random  variables  that  are 
correlated  from  one  time  step  to  another.  When  MIX  = 0 the  process  is 
pure  Gaussian,  and  when  MIX  = 1 the  process  is  pure  lambda-sigma  jump. 

The  example  calculations  use  MIX  = 0.5  which  causes  a mixture  of  the 
two  random  processes. 

SDV(KR)  is  the  standard  deviation  of  fluctuation  phenomenon  KR. 
This  parameter  is  the  primary  method  for  introducing  "modeling  noise"  into 
the  simulation.  The  classification  performance  can  be  degraded  by 
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increasing  SDV.  A value  of  SDV  = 3 dB  is  used  for  four  of  the  five 
random  processes;  SDV  = 1 dB  is  used  for  the  modulation  levels. 

TAU ( KR ) is  the  relaxation  time  of  the  Gaussian  process,  and 
TAU  is  also  the  mean  time  between  jumps  in  the  lambda-sigma  jump  process. 
The  example  uses  TAU  = 3 min  for  all  five  random  processes. 

7 . User-Defined  Function 

FLL( J.KT. ST. DT.ASP)  is  a function  that  calculates  the  Lofar  and 
Demon  line  levels.  As  currently  programmed,  it  is  a dummy  function  that 
sets  FLL  equal  to  the  input  parameter  PLL(J,KT).  A more  sophisticated 
routine  would  include  an  empirically  derived  equation  that  functionally 
relates  line  level  to  target  speed  ST,  depth  DT,  and  aspect  angle  ASP. 

FBB(KT, FRQ, ST, DT, ASP)  is  a function  that  calculates  the  broad- 
band source  spectrum  level  for  target  Type  KT  at  center  frequency  FRQ. 

As  currently  programmed,  the  function  simply  interpolates  the  level/ 
frequency  table,  PBB(NB,KT)/FQB(NB,KT),  to  derive  the  spectrum  level  at 
an  arbitrary  center  frequency.  A more  sophisticated  routine  would  make 
the  spectrum  level  a function  of  target  speed  ST,  depth  DT,  and  aspect 
angle  ASP,  in  addition  to  frequency  and  target  type. 

FNN ( LA . FRQ , SO , DO , ANG ) is  a function  that  calculates  the  broadband 
noise  spectrum  level  at  a center  frequency,  FRQ,  at  the  output  of  Array 
LA.  The  function  first  interpolates  the  level/ frequency  table,  PNN(NF,LA)/ 
FQN(NF,LA),  to  derive  the  spectrum  level  outside  the  array.  Then  the 
directivity  index  is  calculated  for  either  circle  or  line  arrays  by  using 
the  array  dimensions  DH(LA)  and  DV(LA),  and  numbers  of  hydrophones  HN(LA) 
and  VN(LA).  (See  the  source  listing  of  Function  FNN  in  Appendix  A for 
the  exact  method.)  The  noise  level  at  the  array  output  is  then  computed 
by  subtracting  the  directivity  index  from  the  interpolated  noise  level. 

A more  sophisticated  routine  would  make  the  output  noise  level  a function 
of  array  speed  SO,  array  depth  DO,  and  the  point  angle  of  the  sonar 
beam  ANG. 


23 


— - • ^ ^ -ass ^ 


£ 


FAZ(FRQ, RNG,  IX), DT)  is  a function  that  calculates  the  propagation 
loss  at  a center  frequency,  FRQ,  and  range,  RNG.  The  currently  programmed 
function  has  a spreading  term: 

66  + 17  log(RNG) 

and  a frequency  attenuation  term: 

0.08  * RNG  * (FRQ/1,000)L’4  . 

A more  sophisticated  routine  would  contain  a table  look-up  that  included 
the  effects  of  array  depth,  DO,  and  target  depth,  DT,  in  addition  to 
frequency  and  range.  In  the  description  of  the  MUSAC  II  methodology,  two 
propagation  loss  functions  are  possible:  one  a simulation  of  the  real 
propagation  loss,  the  other  a simulation  of  the  sonarman's  expectation  of 
the  propagation  loss.  The  real  function  would  be  used  to  calculate  the 
hypothesized  signal  level  statistics.  The  computer  model  was  simplified 
by  using  FAZ  for  both  the  real  and  expected  propagation  loss  functions. 

FBM(J, FRQ, BRG.ANG)  is  a function  that  calculates  the  beam 
pattern  ratio  for  Array  LA(J)  at  center  frequency  FRQ  when  the  target  is 
on  bearing  BRG  and  the  beam  is  pointed  at  angle  ANG  (both  angles  are 
measured  from  the  course  vector  associated  with  the  array).  As  a way  of 
limiting  the  computations,  the  pointing  angle,  ANG,  is  restricted  to 
equal  the  target  bearings.  Thus  if  there  are  two  target  tracks,  Sub- 
routine FBM  is  used  four  times  each  time  step:  the  beam  is  pointed  at 
Track  1 and  the  response  from  Tracks  1 and  2 is  calculated;  the  beam 
is  then  pointed  at  Track  2 and  the  response  from  Tracks  1 and  2 is  again 
calculated.  The  FBM  routine  contains  four  algorithms,  one  for  each 
combination  of  two  types  of  arrays  (circle,  line)  and  two  types  of  signals 
(narrowband,  broadband).  The  beamwidth  of  a circle  array  remains  constant 
over  pointing  angle,  where  as  the  beamwidth  of  a line  array  increases  when 
the  beam  pointing  angle  approaches  the  direction  angle  of  the  line  array 
(endfire).  The  narrowband  algorithm  calculates  the  nulls  in  the  beam 
pattern.  The  broadband  algorithm  uses  a flat  response;  it  simulates  the 


process  of  averaging  over  frequency,  a process  that  blends  the  single  - 
frequency  beam  pattern  structure  into  a smooth  response  function  over 

2 

angle.  The  basic  beam  pattern  in  all  four  cases  is  a simple  (sin(x)/x) 
response  function.  Details  of  the  function  are  found  in  the  source 
listing  of  Appendix  A. 

FSL(LA, FRQ)  is  a function  that  calculates  the  beam  pattern  ratio 
for  the  reference  channel  of  Array  LA  at  center  frequency  FRQ.  The  value 
of  the  function  is  usually  the  sidelobe  ratio,  SL(LA);  however,  if  the 
mainlobe  of  an  array  is  very  large,  then  the  reference  ratio  may  be  larger 
than  the  sidelobe  input  parameter.  The  reference  beam  pattern  is  used  in 
calculating  the  statistics  of  the  reference  channel.  Because  target 
signals  are  included  in  the  reference  calculations,  the  model  properly 
simulates  sidelobe  masking  of  signals  in  the  mainbeam. 

B.  Output 

The  results  of  the  MUSAC  II  model  are  (1)  range  and  bearing  lists, 

(2)  average  values  of  posterior  probabilities,  and  (3)  probabilities  of 
preference,  classification,  and  tactical  decision.  Probabilities  are 
estimated  by  computing  percentage  of  Monte  Carlo  replications. 

1 . Range  and  Bearing 

The  first  page  of  the  output  lists  the  range,  RNG(LA,M,N),  and 
relative  bearing,  BRG(LA,M, N),  from  each  array,  LA,  to  each  target  track, 

M,  for  each  time  step,  N.  As  shown  in  Appendix  B,  the  first  column  is 
the  range  from  Array  1 (hull-mounted  array)  to  Track  1;  and  the  second 
column  is  the  range  from  Array  2 (towed  away)  to  Track  1.  Since  the 
single-target  example  uses  only  one  track,  the  next  two  columns  do  not 
contain  values;  these  columns  are  used  if  there  is  a second  target  track. 

The  rows  of  the  range  and  bearing  list  represent  time  steps;  the  first 
row  is  at  time  10  min,  the  second  row  is  at  time  20  min,  and  so  on. 

2.  Average  Posterior  Probability 

At  each  time  step,  the  model  calculates  the  posterior  probability, 
POST(I,N),  that  Hypothesis  I is  true,  given  that  the  data  have  accumulated 
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through  time  step  N.  These  probabilities  are  saved  and  averaged  over  the 
replications;  the  average  posterior  probability,  EPOST(I,N),  is  listed 
in  columns  that  represent  the  hypotheses,  and  rows  that  represent  the 
time  steps.  In  Appendix  B,  there  are  three  columns  for  the  hypotheses 
of  the  single-target  example  (Type  1,  Type  2,  and  Type  3 target)  and  10 
rows  for  the  time  steps  from  10  min  to  100  min.  Figure  2 is  a graph  of 
the  average  posterior  as  a function  of  time. 


The  calculation  of  posterior  probability  includes  effects  of 
(1)  the  true  probabilities  of  feature  detection,  (2)  the  hypothesized 
probabilities  of  feature  detection,  (3)  multiple  tracks,  (4)  multiple 
features  on  the  same  and  different  arrays,  (5)  memory  over  time,  and 
(6)  a priori  probabilities.  The  MUSAC  II  methodology  description  must 
be  consulted  to  understand  how  these  effects  are  combined.  The  computer 
implementation  is  similar  to  the  referenced  methodology  except  that 
(1)  sums  of  logarithms  of  likelihoods  were  used  instead  of  products  of 
likelihoods,  and  (2)  the  total  likelihood  was  scaled  from  1 to  1,000-- 
thus  no  hypothesis  could  be  more  than  1,000  times  more  likely  than 
another.  These  computer  implementation  steps  were  taken  to  avoid  certain 
numeric  problems  in  the  calculation  of  posterior  probabilities. 

The  average  posterior  probability  is  included  as  a simulation 
result  because  the  values  can  give  an  indication  of  feature  detection 
status.  Average  posteriors  that  are  nearly  equal  arise  because  dis- 
tinguishing features  are  not  detected;  dissimilar  average  posteriors 
indicate  that  combinations  of  distinguishing  features  were  detected. 
Posterior  probabilities  are  not  probabilities  of  classification,  and  are 
only  indicative  of  the  classification  results. 

3 . Probability  of  Preference 

The  probability  of  preference,  P0P(I,N),  is  the  percent  of 
replications  for  which  Hypothesis  I was  the  preferred  answer  at  Time 
Step  N.  If  a classification  decision  must  be  made  at  Time  Step  N (not 
before  N,  not  after  N),  then  the  probability  of  preference  would  equal 
the  probability  of  classification.  But  the  model  allows  for  only  one 
classification  decision,  which  can  occur  at  any  time  step;  therefore, 
a different  name  was  coined  for  the  concept  of  step-to-step  classifica- 
tion probabilities.  In  determining  the  probability  of  preference,  the 
I-th  hypothesis  is  preferred  if  P0ST(I,N)  is  larger  than  all  the  other 
posterior  probabilities  at  time  step  N. 
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The  example  output  in  Appendix  B lists  the  probabilities  of  | 

preference  in  columns  that  represent  hypotheses,  and  rows  that  represent 
time  steps  (the  values  are  in  increments  of  0.05  because  20  replications 
were  used).  Figure  3 shows  the  example  probabilities  of  preference  as 
a function  of  time. 


4. 


Probability  of  Classification 

The  probability  of  classification,  PLC(I),  is  the  percent  of 
replications  for  which  the  I-th  hypothesis  was  was  chosen.  Since  the 
probability  of  classification  does  not  contain  any  information  about  when 
the  decision  was  made,  the  time  distribution  for  each  decision,  DISTR(I,N), 
is  also  calculated.  Under  the  condition  that  the  decision  was  ultimately 
Hypothesis  I,  DISTR(I,N)  is  the  probability  that  the  decision  was  made  at 
or  before  Time  Step  N.  The  probabilities  of  classification  and  the  time 
distribution  are  on  page  3 of  the  output.  In  the  single  target  example, 

55  percent  were  correctly  classified,  15  percent  were  classified  as  Type 
2 target,  and  30  percent  were  classified  as  Type  3 target.  Figure  4 shows 
how  the  55  percent  correct  classification  decisions  were  distributed  in 
time . 

The  first  three  pages  of  the  output  in  Appendix  B show  the 
results  of  Run  1,  and  the  next  six  pages  are  the  results  of  two  additional 
runs  of  the  simulation.  The  three  runs  differ  in  the  value  of  the 
parameter  IR  which  is  the  true  hypothesis  index: 

Run  Real  Target 

1 Type  1 

2 Type  2 

3 Type  3 

When  the  results  of  the  three  runs  are  combined,  a matrix  of  the  classi- 
fication probabilities  (in  units  of  percent)  can  be  constructed: 


Real 

Target 

Classified  as 

1 2 3 

1 

55 

15 

30 

2 

10 

65 

25 

3 

5 

10 

85  . 
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CONDITIONAL  PROBABILITY  DISTRIBUTION 


FIGURE  4 DISTRIBUTION  OF  A CLASSIFICATION  DECISION 
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For  example,  when  the  real  target  was  Type  3,  5 percent  of  the  replica- 
tions (1  out  of  20)  were  classified  as  Type  1.  In  the  remainder  of  this 
report,  the  above  matrix  is  called  the  "classification  matrix,"  and  only 
the  essential  information  is  displayed: 


55  15 

10  65 


15  30 

65  25 

10  85 


5 . Probability  of  Tactical  Decision 

The  probability  of  tactical  decision,  PKD(KD),  is  the  percent 
of  replications  for  which  the  KD-th  tactical  alternative  was  chosen. 
Tactical  alternatives  are  related  to  hypotheses  through  the  value  matrix, 
VAL(I,KD).  For  the  output  shown  in  Appendix  B,  the  value  matrix  was  a 
3-by-3  identity  matrix,  and  thus  the  tactical  alternatives  were  the  same 
as  the  hypotheses.  This  is  the  reason  that  the  probability  of  tactical 
decision  output  is  identical  to  the  probability  of  classification  output. 

To  see  the  effect  of  a different  value  structure,  the  value 
matrix  was  changed  to: 


Hypothesis 

Tactical 

Alternative 

1 2 

1 

300 

-500 

2 

-100 

100 

3 

-500 

100  , 

By  making  three  runs  of  the  simulation,  a "tactical  alternative  matrix" 
was  computed;  it  shows  the  percent  of  replications  for  which  Alternative 
1 or  2 was  chosen  for  each  of  three  real-target  conditions: 


Real  Tactical 
Target  Alternative 
1 2 


0 100 
31 


- 


1 

2 

3 


The  tactical  decision  is  based  on  maximizing  engagement  value, 
determined  from  the  value  matrix  and  the  posterior  probabilities.  A 
classification  decision  on  the  best  hypothesis  is  also  made  by  choosing 
the  highest  posterior  probability  at  the  time  of  the  tactical  decision. 
Therefore,  a classification  matrix  was  produced  using  the  new  value 
structure,  and  this  matrix  is  compared  to  the  classification  matrix  of  the 
original  simulation: 


New  Value  Base 

Structure  Case 


50 

15 

35 

55 

15 

30 

0 

70 

30 

10 

65 

25 

0 

0 

100 

5 

10 

85 

Changing  the  value  matrix  from  an  identity  matrix  to  a matrix  representing 
tactical  tradeoffs  changes  the  classification  probabilities,  sometimes  for 
the  better  and  sometimes  for  the  worse. 

The  timing  of  the  tactical  decision  is  also  a simulation  output. 
Figure  5 shows  the  time  distribution  of  tactical  decisions  for  the  case 
where  the  real  target  was  Type  1.  Both  decision  alternatives  are  shown 
on  the  figure;  Alternative  1 is  defined  as  an  "attack"  decision  and  Alter- 
native 2 is  a "no-attack"  decision.  For  a Type  1 target,  the  attack 
decisions  were  more  likely  than  the  no  attack  decisions  (55  percent 
versus  45  percent)  and  they  occurred  at  shorter  ranges  (longer  times)  than 
the  no-attack  decisions. 

6 . Engagement  Measure -of -Effectiveness 

The  final  output  statistic  is  a single  number  called  the 
engagement  measure-of -effectiveness  (MOE).  It  is  the  average  value  of  the 
final  tactical  decision  based  on  the  true  hypothesis;  it  is  averaged 
over  all  the  replications: 

MOE  = [Z  VAL( IR,KDF)] /NREP  , 
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where  IR  is  the  true  hypothesis,  KDF  is  the  final  tactical  decision, 
and  NREP  is  the  number  of  replications.  When  the  value  matrix  is  the 
identity  matrix,  the  MOE  is  just  the  probability  of  correct  classification 


NO-ATTACK 

DECISION 

(45%) 


■ ATTACK 
DECISION 
(55%) 


REAL  TARGET:  TYPE  1 


TIME  — minutes 


FIGURE  5 DISTRIBUTION  OF  TACTICAL  DECISIONS 


With  the  altered  value  structure  of  the  previous  section,  the  engagement 
MOEs  for  the  three  runs  were: 

Real 

Target  MOE 

1 -60 

2 100 

3 100 

Even  though  the  correct  classification  probabilities  for  Type  2 and 
Type  3 targets  were  different  (70  percent  versus  100  percent)  the  MOEs 
were  the  same,  and  thus  engagements  with  either  of  these  targets  types 
is  equally  "valuable". 


C.  Variations  on  the  Base  Case 

The  following  sections  describe  the  sensitivity  of  single-target 
simulation  results  when  selected  parameters  are  varied  from  their  base- 
case  values. 

1 . Lofar  and  Demon 

The  model  was  run  with  only  Lofar  features  and  then  run  again 
with  only  Demon  features.  These  two  cases  were  then  compared  to  the 
base  case  classification  matrix  in  which  both  Lofar  and  Demon  were  used: 


Lofar  Only 


Demon  Only 


Both  Lofar 
and  Demon 


70 

10 

20 

65 

20 

15 

55 

15 

30 

15 

85 

0 

25 

55 

20 

10 

65 

25 

10 

25 

65 

15 

15 

70 

5 

10 

85 

When  the  real  target 

was  Type 

1 

or  2, 

the 

use  of 

both 

Lofar  and  Demon 

features  degraded  the  correct  classification  probabilities  relative  to 
the  Lofar-only  or  Demon-only  cases.  However,  when  the  target  was  Type  3, 
the  use  of  both  sets  of  features  increased  the  probability  of  correct 
classification.  Adding  more  features  to  aid  in  classifying  a target  does 
not  necessarily  improve  the  classification  performance;  in  fact,  adding 
more  features  may  degrade  the  performance. 
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2. 


Decision  Threshold 


When  the  decision  threshold,  FVAL,  was  varied,  the  classifica- 
tion matrices  changed  and  the  time  distribution  of  classification 
decisions  also  changed.  For  example,  the  classification  matrices  were 
computed  for  three  valuta  of  the  threshold: 


FVAL  = 

0.8 

FVAL 

= 0. 

95 

FVAL 

= 0. 

99 

35 

25 

40 

55 

15 

30 

50 

15 

35 

10 

60 

30 

10 

65 

25 

0 

80 

20 

5 

20 

75 

5 

10 

80 

10 

10 

80 

In  general,  by  reducing  the  threshold,  there  are  fewer  correct 
classification  decisions  but  they  are  made  sooner.  The  time  distribution 
of  the  correct  classification  decisions  for  the  three  FVAL  runs  is  shown 
in  Figure  6. 

3 . Standard  Deviation 

The  random  process  standard  deviation  vector,  SDV(KR),  was 
varied  from  the  base  case  value  by  subtracting  and  adding  1 dB: 

Variation  1:  SDV  = 2,  2,  0,  2,  2 

Base  case:  SDV  = 3,  3,  1,  3,  3 

Variation  2:  SDV  = 4,  4,  2,  4,  4 

The  SDV  components  are  Lofar,  BBand,  Demon,  PLoss,  and  Noise,  respectively. 

The  resulting  classification  matrices  were: 


Variation 

1 

Base 

Case 

Variation 

2 

60 

15 

25 

55 

15 

30 

50 

20 

30 

5 

70 

25 

10 

65 

25 

15 

65 

20 

10 

5 

85 

5 

10 

85 

15 

20 

65 

When  SDV  became  larger,  more  classification  mistakes  were  made.  Within 
limits,  the  SDV  parameter  may  be  used  to  adjust  the  model  results  to 
correspond  with  experimental  data  points.  Then  the  model  can  predict 
classification  results  for  cases  not  covered  in  the  experiment. 


35 


Ill  DOUBLE-TARGET  EXAMPLE 


This  chapter  describes  how  MUSAC  II  can  simulate  a scenario  in 
which  there  are  two  targets.  Many  of  the  input  parameters  used  in  the 
single-target  example  are  used  again  in  the  double-target  example;  only 
the  changes  to  the  single-target  example  are  discussed.  Appendix  C con- 
tains: (1)  a listing  of  the  double-target  base-case  input  parameters 

and  (2)  a listing  of  the  base-case  output. 

The  double-target  scenario  involves  the  passive  acoustic  classifi- 
cation of  two  surface  ships  by  a submarine  using  Demon  information  from 
a hull-mounted  array  and  Lofar  information  from  a towed  array.  The  sce- 
nario geometry  is  shown  in  Figure  7.  Initial  range  separation  between 
the  submarine  and  the  two  target  ships  is  about  40  nmi.  The  surface 
ships  are  5 nmi  apart  and  travel  on  a near-collision  course  with  the 
submarine.  During  the  1.5-hour  engagement  the  submarine  must  classify 
the  targets  as  one  of  nine  possible  target  configurations.  Figure  8 
shows  the  relative  bearing  to  the  two  target  tracks  as  a function  of 
time.  The  initial  target  separation  is  about  6 degrees  and  the  final 
separation  about  36  degrees. 

A.  Base  Case 

The  double-target  base  case  example  is  used  to  demonstrate  the 
multitarget  capability  of  MUSAC  II  and  to  provide  results  for  comparison 
with  cases  in  which  selected  input  parameters  are  varied. 

1 . Input 

Only  five  input  parameters  need  to  be  changed  to  turn  the 
single-target  example  into  a double-target  example. 

IMAX  is  increased  from  three  hypotheses  to  nine  hypotheses, 


which  are  detailed  below. 


NMAX  is  increased  from  1 track  to  2 tracks;  in  this  way,  two 
targets  can  be  simulated. 

KT(I ,M)  is  increased  from  a 3-by-l  array  to  a 9-by-2  array. 
The  nine  hypotheses  are  defined  by  the  KT  array  as  follows: 


Target 

Type 

: on  Track 

Hypothesis 

1 

2 

1 

1 

1 

2 

1 

2 

3 

1 

3 

4 

2 

1 

5 

2 

2 

6 

2 

3 

7 

3 

1 

8 

3 

2 

9 

3 

3 

For  example.  Hypothesis  8 states  that:  "Target  Type  3 is  on  Track  1 and 
Target  Type  2 is  on  Track  2."  For  small  numbers  of  target  types  and 
tracks,  the  above  combinational  method  of  constructing  hypotheses  can  be 
used.  However,  when  the  scenario  is  complex,  the  analyst  must  reduce 
the  number  of  hypotheses  by  excluding  the  ones  with  very  low  a priori 
probability. 

IR  is  the  hypothesis  number  that  represents  the  real  target 
configuration.  Hypothesis  3 was  chosen  as  reality  for  the  base  case 
(IR  = 3).  Therefore,  a target  of  Type  1 is  on  Track  1 and  a target  of 
Type  3 is  on  Track  2. 

KDMAX  is  changed  from  a total  of  3 to  9 tactical  alternatives 
so  that  the  tactical  alternatives  remain  the  same  as  the  nine  classifi- 
cation alternatives.  As  a variation  on  the  base  case,  an  engagement 
value  structure  that  is  not  an  identity  matrix  is  used. 


2 . Output 


Appendix  C shows  the  results  of  the  double- target  base  case. 
The  range  and  relative  bearing  versus  time  are  listed  for  four  columns 
of  array/ track  combinations. 


Column 


Array 


1 

2 

3 

4 


1 

2 

1 

2 


Track 


1 

1 

2 

2 


The  relative  target  bearing  from  Array  1 to  Tracks  1 and  2 (Columns  1 
and  3)  were  previously  shown  in  Figure  8. 

The  average  posterior  probability  and  the  probability  of  a 
preference  are  listed  in  9 columns  and  10  rows,  which  relate  to  the  9 
hypotheses  and  10  time  steps.  Figure  9 shows  the  probability  of  prefer- 
ence for  Hypothesis  3 (the  correct  hypothesis)  plotted  as  a function  of 
time  on  the  lowest  curve.  Two  other  curves  are  shown  for  comparison. 

The  middle  curve  is  the  sum  of  the  preference  probabilities  for  Hypothe- 
ses 1,  2,  and  3;  it  is  the  probability  that  target  Type  1 is  on  Track  1 
and  that  any  target  type  was  on  Track  2.  The  top  curve  is  the  sum  of  the 
preference  probabilities  for  Hypotheses  1,  2,  3,  4,  and  7;  it  is  the 
probability  that  target  Type  1 was  present.  Figure  9 demonstrates  that 
probabilities  may  be  added  together  to  construct  higher-order  hypotheses 
using  the  elemental  hypotheses  of  the  model. 

The  probability  of  classification  (in  percent)  for  the  double- 
target base  case  was  computed  as: 


Real 

Conf iguration 

Classif ied 

1,1  1,2  1,3  2,1 

Conf iguration 

2,2  2,3  3,1  3,2  3,3 

1,3 

35  5 20  5 

0 5 10  0 20 

Correct  classification  occurred  for  only  20  percent  of  the  replications, 
and  classifying  the  correct  target  type  on  Track  1 (1,1  1,2  1,3)  occurred 


for  60  percent  of  the  replications.  A classification  matrix  could  have 
been  computed  by  performing  nine  runs  and  changing  the  true  configuration 
each  run.  The  above  result  is  one  row  of  the  classification  matrix. 

Only  two  double-target  classification  decisions  were  made  be- 
fore the  last  time  step,  whereas  all  20  single-target  decisions  were  made 
before  the  last  time  step.  Classifying  multiple  target  configurations  is 
harder  than  single  targets  when  the  same  decision  threshold  is  used. 


(-500)  + (100)  = (-400).  Of  course,  other  methods  can  be  devised  to 
produce  a value  matrix  that  reflects  the  value  of  a tactical  action  taken 
against  a multitarget  group. 

The  computed  probabilities  of  classification  (in  percent)  for 
the  two  cases  were: 


Classified 

Conf iguration 

Base 

Case 

New 

Value 

Structure 

1,1 

35 

40 

1,2 

5 

5 

1,3  (true) 

20 

35 

2,1 

5 

0 

2,2 

0 

0 

2,3 

5 

0 

3,1 

10 

10 

3,2 

0 

0 

3,3 

20 

10 

The  result  of  using  the  new  value  structure  was  an  increase  in  correct 
classifications . 

The  probability  of  tactical  decisions  was  65  percent  for  the 
attack  alternative  and  35  percent  for  the  no-attack  alternative.  The 
median  time  for  an  attack  decision  was  about  40  min  and  the  median  time 
for  a no-attack  decision  was  about  15  min.  The  time  distribution  of 
decision-making  under  the  new  value  structure  was  sharply  different  from 
the  distribution  for  the  base  case,  where  almost  all  decisions  were  de- 
layed until  the  end  of  the  engagement. 

2 . Decision  Threshold 

Figure  10  shows  how  the  engagement  MOE  changed  as  a function 
of  the  decision  threshold  when  the  tactical  value  matrix  of  the  last  sec- 
tion was  used.  The  maximum  engagement  value  is  attained  when  the  deci- 
sion is  to  attack,  and  the  minimum  when  the  decision  is  not  to  attack. 
This  graph  of  engagement  MOE  versus  a system  parameter  demonstrates  how 
system  analyses  may  be  performed  with  the  MUSAC  II  simulation  model. 
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DECISION  THRESHOLD  — FVAL 

FIGURE  10  MOE  SENSITIVITY  TO  DECISION  THRESHOLDS 
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PROGRAM  MUSAC2( OUTPUT, TAPE6=OUTPUT ) 


C 

C 

c 


c 

INPUT 

c 

XX 

GV 

c 

X* 

NREP 

c 

X X 

NMAX 

c 

XX 

TS 

c 

X X 

I MAX 

c 

XX 

MMAX 

c 

XX 

KT ( I , M ) 

c 

X X 

IR 

c 

X X 

PR1 OR( I ) 

c 

X X 

KDMAX 

c 

XX 

VAL ( 1 , KD) 

c 

X X 

FVAL 

c 

XX 

XT  ( M ) 

c 

XX 

YT  ( M ) 

c 

X X 

CT ( N, M) 

c 

X X 

ST ( N , M ) 

c 

X X 

DT( N, M) 

c 

X X 

JMAX 

c 

X X 

NF  ( J ) 

c 

X X 

KF  ( J ) 

c 

XX 

FRQ( J ) 

c 

X X 

PLL ( J , KT ) 

c 

X X 

PLL( J , KT ) 

c 

X X 

PBB( NB , KT ) 

c 

X X 

FQB( NB, KT ) 

c 

XX 

LAMAX 

c 

XX 

LA  ( J ) 

c 

X X 

KA(LA) 

c 

XX 

DH(LA) 

c 

X X 

DV ( LA ) 

c 

X X 

HN(LA) 

c 

X X 

VN(LA) 

c 

X X 

SL(LA) 

c 

X X 

XO ( LA ) 

c 

X X 

YO(LA) 

c 

X X 

CO( N, LA) 

c 

X X 

SO(N, LA) 

c 

X X 

D0( N, LA) 

c 

XX 

NA(N, LA) 

c 

X X 

PNN(NF, LA) 

c 

X X 

FQN( NF, LA) 

c 

X X 

LP  ( J ) 

c 

X X 

WTH(LP) 

c 

X X 

SCR ( LP ) 

c 

X X 

FCS(LP) 

c 

X X 

TOT ( LP) 

c 

X X 

DET(LP) 

c 

X X 

MIX(KR) 

c 

X X 

SDV(KR) 

c 

X X 

TAU(KR) 

RANDOM  NUMBER  GENERATIVE  VALUE 
NUMBER  OF  REPLICATIONS 
NUMBER  OF  TIME  STEPS 
TIME  STEP  DURATION  (MIN) 

NUMBER  OF  HYPOTHESES 

NUMBER  OF  TRACKS 

HYPOTHESIZED  TARGET  TYPE 

TRUE  HYPOTHESIS  NUMBER 

PRIOR  PROBABILITY  OF  HYPOTHESIS 

NUMBER  OF  TACTICAL  DECISION  ALTERNATIVES 

VALUE  OF  DECISION 

DECISION  THRESHOLD  OF  RATIO  ( MAX . EXP . VALUE )/( EXP . MAX . VALUE ) 
TARGET  INITIAL  POSITION  ( NMI ) 

TARGET  INITIAL  POSITION  (NMI) 

TARGET  COURSE  (DEG) 

TARGET  SPEED  (KT) 

TARGET  DEPTH  (FT) 

NUMBER  OF  FEATURES 
FEATURE  OFF/ON  (0=0FF  1=ON> 

FEATURE  TYPE  ( 1 =LOFAR  2=BBAND  3=DEMON) 

CENTER  FREQUENCY  (HZ) 

LOFAR  LINE  LEVEL  ( DB  UPA2  1 YD ) ALSO 
DEMON  MODULATION  LEVEL  (DB) 

BROADBAND  SOURCE  SPECTRUM  OF  TARGET  ( DB  UPA2/HZ  1 YD ) 
FREQUENCY  POINTS  FOR  PBB  (HZ) 

NUMBER  OF  ARRAYS 
ARRAY  NUMBER 

ARRAY  TYPE  (1=CIRCLE  2=LINE) 

HORIZONTAL  ARRAY  SIZE  (FT) 

VERTICAL  ARRAY  SIZE  (FT) 

HORIZONTAL  NUMBER  OF  HYDROPHONES 
VERTICAL  NUMBER  OF  HYDROPHONES 
SI  DELOBE  RATIO 
ARRAY  INITIAL  POSITION  (NMI) 

ARRAY  INITIAL  POSITION  (NMI) 

ARRAY  COURSE  (DEG) 

ARRAY  SPEED  (KT) 

ARRAY  DEPTH  (FT) 

ARRAY  OFF/ON  (0=0FF  1=ON) 

BROADBAND  NOISE  OUTSIDE  ARRAY  ( DB  UPA2/HZ ) 

FREQUENCY  POINTS  FOR  PNN  (HZ) 

PROCESSOR  NUMBER 
BANDWIDTH  (HZ) 

SCAN  RATE  ( NUMBER/MI N) 

NUMBER  OF  FEATURE  CELLS  PER  SCAN 
TOTAL  TIME  OF  FEATURE  INTEGRATION  (MIN) 

DETECTION  THRESHOLD  (NUMBER  OF  SIGMAS) 

MIXING  CONSTANT  (0.=GAUSS  1.=JUMP)  KR  1 = LOFAR  4=NOISE 
STD  DEV  OF  RANDOM  PROCESS  ( DB ) 2=BBAND  S=PLOSS 

RELAXATION  TIME  OF  RANDOM  PROCESS  (MIN)  3=DEMON 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


**  FLL (J,KT,ST,DT, ASP ) 

**  FBB(KT, FRQ, ST, DT, ASP) 
**  FNN ( LA , FRQ , SO , DO , ANS ) 
**  FAZ(FRQ, RNG, DO, DT) 

**  FBM( J, FRQ, BRG, ANG) 

**  FSL ( LA, FRQ ) 


LOFAR  AND  DEMON  SOURCE  LEVEL 
BROADBAND  SOURCE  SPECTRUM  LEVEL 
NOISE  SPECTRUM  LEVEL  AT  ARRAY  OUTPUT 
REAL  PROPAGATION  LOSS  <DB> 
BEAMPATTERN  RATIO 
REFERENCE  BEAMPATTERN  RATIO 


c 

INTERNAL 

c 

X K 

I 

HYPOTHESIS  NUMBER 

c 

* * 

J 

FEATURE  NUMBER 

c 

* * 

K 

LOOK  ANGLE  NUMBER 

c 

* X 

M 

TRACK  NUMBER 

c 

X X 

N 

TIME  STEP  NUMBER 

c 

X X 

NN 

REPLICATION  NUMBER 

c 

X X 

KR 

KIND  OF  RAMDOM  PROCESS 

c 

X X 

3 = M0DULAT I ON  . 

c 

XX 

KD 

TACTICAL  DECISION  NUMBI 

C **  TAV(LP) 

C *«  NUM(LP) 

C *«  MEM(LP) 

C ««  RNG( LA, M, N) 

C **  BRG ( LA, M, N) 

C **  ASP ( LA , M, N ) 

C **  ANG( LA, K, N) 

C «»  TLL(J,M,N) 

C **  TBB(  J , M,  II) 

C **  TNN( J , K , N) 

C **  TAZ ( J , M, N ) 

C **  TBM(JKMN) 

C **  TSLIJ) 

C **  XS(IJKN) 

C **  XV( I JKN ) 

C **  XSP(IJKN) 

C **  XDSP(IJKN) 

C **  JUMP (KR, J , M) 

C **  GAUSS (KR, J , M) 

C **  BIAS(KR) 

C **  PDETZ ( J , K ) 

C *«  PDET ( I , J, K ) 

C **  LIKE( I , N) 

C **  POST ( I , N ) 

C **  I B ( N ) 

C **  VMAX ( I ) 

C **  KDF 
C **  NSTOP 
C INTERNAL  FUNCTIONS 
C **  FLBC ( N, K ) 

C *«  PROBCMU, SIG, THRESH) 

C **  RNORM( 0 . ) 

C **  RANDOM( MI  X , SDV , TAU, TS , 
C «*  JUMP,  GAUSS, DELTA) 


1=LINE  2=BROADBAND 


AVERAGING  TIME  PER  FEATURE  CELL  (SEC) 

NUMBER  OF  SCANS  PER  TIME  STEP 

NUMBER  OF  TIME  STEPS  IN  LIKELIHOOD  CALCULATION 
RANGE  ( NMI ) 

RELATIVE  BEARING  (DEG) 

TARGET  ASPECT  (DEG) 

SONAR  LOOK  ANGLE  (DEG) 

LINE  LEVEL  TABLE  (DB) 

BROADBAND  SPECTRUM  TABLE  (DB) 

NOISE  TABLE  (DB) 

REAL  PROP  LOSS  TABLE  (DB) 

BEAM  PATTERN  TABLE  (RATIO) 

SIDE  LOBE  TABLE  (RATIO) 

HYPOTHETICAL  LOFAR  SIGNAL 
HYPOTHETICAL  LOFAR  SIGNAL  SQUARED 
HYPOTHETICAL  BROADBAND  SIGNAL 
HYPOTHETICAL  DEMON  SIGNAL 

LAST  VALUE  OF  JUMP  PROCESS  RANDOM  VARIABLE 
LAST  VALUE  OF  GAUSSIAN  PROCESS  RANDOM  VARIABLE 
BIAS  ADDED  TO  COMPENSATE  FOR  RANDOM  PROCESS  (DB) 
PROBABILITY  OF  DETECTION  OF  TRUE  DATA 
PROBABILITY  OF  DETECTION  OF  HYPOTHETICAL  DATA 
LIKELIHOOD  OF  DATA 

POSTERIOR  PROBABILITY  OF  HYPOTHESIS 
MAXIMUM  POSTERIOR  HYPOTHESIS 
MAXIMUM  OF  VAL( I , KD)  OVER  KD 
FINAL  TACTICAL  DECISION 
TIME  STEP  OF  FINAL  DECISION 


LOG  OF  BINOMIAL  COEFFICIENT 
PROB  NORMAL  RAN. VAR.  GE  THRESHOLD 
RANDOM  NORMAL  DEVIATE  (MU=0  SIG=1) 
CALCULATES  ZERO -MEAN  RANDOM 

DEVIATE  FROM  MIXED  PROCESS 
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r i 


1 

OUTPUT 

**  RNG(LA,M,N) 

**  BRG( LA,M,N) 

**  EPOST ( I , N ) 

*#  POP( I , N) 

**  PCL( I ) 

**  D I STR ( 1 , N ) 

**  PKD(KD) 

**  H I STO ( KD, N ) 

**  VMOE 


RANGE  (NM1  > 

RELATIVE  BEARING  (DEG) 

AVERAGE  POSTERIOR  PROBABILITY 
PROBABILITY  OF  PREFERENCE 
PROBABILITY  OF  CLASSIFICATION 
DISTRIBUTION  OF  CLASSIFICATION  DECISION 
PROBABILITY  OF  TACTICAL  DECISION 
DISTRIBUTION  OF  TACTICAL  DECISION 
ENGAGMENT  MEASURE  OF  EFFECTIVENESS 


COMMON  /A/  GV, NREP, NMAX 

COMMON  /B/  I MAX, MMAX, JMAX 

COMMON  /C/  LAM  AX,  KDMAX,  I R 

COMMON  /D/  LP ( 1 2 ) , LA ( 1 2 ) , NA ( 20 , 2 ) , KA ( 2 ) 

COMMON  /E/  KF( 12) , NF( 12) ,KT( 10, 2)  , TS 

COMMON  /F/  X0( 2 ) , YO ( 2 ) , C0( 20,2), SO ( 20, 2 ) , DOC  20, 2 ) 

COMMON  /G/  XT { 2 ) , YT ( 2 ) , CT( 20, 2 ) , ST ( 20,  2 ) , DT ( 20,  2 ) 

COMMON  /H/  WTH ( 5 ) , TAV ( 5 ) , NUM ( 5 ) , DET ( 5 ) , MEM ( 5 ) 

COMMON  /I/  DHC  2) , DV ( 2) , HN( 2) , VN( 2) , SL ( 2) 

COMMON  /J/  FRQ<  12)  , PLU12,  3) 

COMMON  /K/  PBB (6,3), FQB (6,3) 

COMMON  /L/  PNN (6,2), FQN( 6,2) 

COMMON  /M/  M I X ( 5 ) , SDV (5),TAU(5),BI AS ( 5 ) 

COMMON  /N/  PRIORt 10) , VAL( 1 0, 10) , FVAL, VMAX( 10) 

COMMON  /O/  RNG ( 2, 2, 20) , BRG( 2, 2, 20 ) , ASP (2,2,20) 

COMMON  /P/  TLL( 12, 2, 20) , TBB( 12, 2, 20) , TNN( 12, 2, 20) 

COMMON  /Q/  TAZ( 1 2, 2, 20) , TBM( 960)  , TSL ( 1 2 ) 

COMMON  /R/  JUMP(5, 12,2), GAUSS(5,  12,2) 

COMMON  /S/  PD£TZ( 12, 2) , PDETC 10,  12,  2) 

COMMON  /T/  POSTC 10, 20) , IB(20) 

COMMON  /U/  NSTOP, KDF 

COMMON  /V/  XS (4800 ) , XV (4800 ) , XSP (4800 ) , XDSP ( 4800 ) 

DIMENSION  EPOST ( 10, 20) , DISTRt 10,20), TOTAL (10), PCL( 10)  , 

+ HI STO ( 10, 20) ,SUM( 10) , PKD( 10), POP( 10,20) 

REAL  MIX 

CALL  INPUT 
DO  90  IR=1 ,3 
CALL  TABLE 
CALL  SAVE 
DO  12  1 = 1,1  MAX 
TOTAL ( I ) =0 . 

DO  12  N= 1 , NMAX 
EPOST ( I , N) =0. 

POP( I , N) =0 . 

12  D I STR ( I , N ) =0 . 

DO  14  KD= 1 , KDMAX 
SUM ( KD ) =0 . 

I 


o o o o o 


DO  14  N=  1 , NMAX 
14  H I STO( KD, N ) =0 . 

VM0E=0. 

**  REPLICATION  NN-LOOP 
CALL  RANSET(GV) 

DO  50  NN= 1 , NREP 
**  INITIALIZE  RANDOM  PROCESSES 
DO  1 0 KR= 1 , 5 
DO  10  J=  1 , JMAX 
DO  10  M= 1 , MMAX 

JUMP ( KR , J,  M ) =RNORM ( 0 . )«SDV(KR) 
10  GAUSS ( KR , J,M) =RNORM ( 0 . )«SDVCKR) 

**  TIME  STEP  N-LOOP 
DO  30  N=1 , NMAX 
CALL  DETECT! N) 

CALL  BAYES! N) 

DO  25  1=1,1  MAX 

25  EPOST! I , N) =EPOST! I , NJ+POST! I , N) 
I BN= I B ( N ) 

POP! I BN, N ) =POP ! I BN , N ) + 1 . 

30  CONTINUE 

CALL  DECIDE 
I I = IB(NSTOP) 

TOTAL! I I ) =TOTAL! I I ) +1  . 

DO  40  N=NST0P, NMAX 
40  DISTR! I I ,N)=DISTR( I I , N ) +1  . 

SUM! KDF) =SUM IKDF) +1  . 

DO  45  N=NSTOP, NMAX 
45  HI STO! KDF, N ) =HI STO (KDF, N) +1 . 

VMOE= VMOE+VAL ! I R, KDF ) 

50  CONTINUE 

DO  60  1=1,  I MAX 
PCL! I )=TOTAL! I 1/NREP 
DO  60  N= 1 , NMAX 
EPOST! I , N) =EPOST! I , N) /NREP 
POP! I , N) =POPl I , N) /NREP 
IF! TOTAL! I ) . EG . 0 . ) GO  TO  60 
DISTR! I , N) =DI STR! I , N) /TOTAL!  I ) 
60  CONTINUE 

DO  70  KD= 1 , KDMAX 
PKDIKD) =SUM (KD) /NREP 
I F( SUM!KD)  . E Q.O.)  GO  TO  70 
DO  65  N= 1 , NMAX 

65  HI  STO (KD, N ) =HI STO(KD, N ) /SUM (KD) 
70  CONTINUE 

VMOE=VMOE/NREP 


WRI 

TE  ( 6, 

100) 

WRI 

T:(6, 

150) 

RNG 

WRI 

TE!  6, 

120) 

WRI 

TE  ( 6 , 

150) 

BRG 

WRI 

TE  ( 6, 

200) 

50 


90 


WRI TE ( 6, 250 ) EPOST 
WR I TE( 6, 270) 
WRITEI6, 230) 

WRITE (6, 300) 

WRI TEI6, 230) 

WR I TE( 6, 400) 

WR I TE ( 6, 250 ) DISTR 
WR I TE( 6, 500) 
WRITEC6, 250) 

WRI TE!6, 600) 

WR I TE! 6 , 250) 

WRI TE(6, 700) 
CONTINUE 


POP 


PCL 


PKD 

HISTO 

VMOE 


100 

120 

200 

270 

300 

400 

500 

600 

700 

150 

250 


FORMAT ( 1 HI , *RANGE* ) 

FORMAT!//*  BEARING*) 

FORMAT! 1H1 , *AVERAGE  POSTERIOR  PROBABILITY*) 
FORMAT!//*  PROBABILITY  OF  PREFERENCE*) 
FORMAT! 1 HI,  * PROBABILITY  OF  CLASSIFICATION*) 


FORMAT!/*  DISTRIBUTION  OF 
FORMAT!//*  PROBABILITY  OF 
FORMAT!/*  DISTRIBUTION  OF 
FORMAT!//*  ENGAGMENT  MOE 
FORMAT ( 5X , 4F 1 0 . 2 ) 

FORMAT ! 5X, 10F10. 2) 

END 


CLASSIFICATION  DECISION*) 
TACTICAL  DECISION*) 
TACTICAL  DECISION*) 

=*, F7. 3) 
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SUBROUTINE  INPUT 
C 


COMMON  /A/  GV, NREP, NMAX 

COMMON  /B/  I MAX , MMAX , JMAX 

COMMON  / C/  LAMAX, KDMAX, I R 

COMMON  /D/  LP ( 1 2 ) , LA  I 1 2 ) , NA (20,2), KA ( 2 ) 

COMMON  /E/  KF( 12) , NFI 12) ,KT( 10, 2) , TS 
COMMON  /F/  X0( 2) , YO ( 2 ) , C0( 20, 2) , SOI  20, 2) , DO ( 20, 2 ) 
COMMON  /G/  XT ( 2) , YT ( 2 ) , CT( 20, 2) , ST I 20, 2) , DT I 20, 2 ) 
COMMON  /H/  WTH( 5 ) , TAV ( 5 ) , NUM( 5 ) , DET( 5 ) , MEM  I 5 ) 

COMMON  /I/  DH ( 2 ) , DV ( 2 ) , HN ( 2 ) , VN ( 2 ) , SL ( 2 ) 

COMMON  /J/  FRQ( 12) , PLL( 12, 3) 

COMMON  /K/  PBB( 6,3), FQB( 6,3) 

COMMON  /L/  PNN (6,2), FQN (6,2) 

COMMON  /M/  M I X ( 5 ) , SDV (5),TAU(5),BI AS ( 3 ) 

COMMON  /N/  PRIORI  10) ,VAL(10, 10) , FVAL, VMAX1 10) 

DIMENSION  SCRC  5) , FCS ( 5 ) , TOT I 5) 

REAL  MIX 

DATA  GV  /808 . / 

DATA  NREP  /20/ 

DATA  NMAX  /10/ 

DATA  TS  /10./ 

DATA  I MAX  /3/ 

DATA  MMAX  /I/ 

DATA  KT  /I , 2, 3,  1 7*0/ 

DATA  I R / 1 / 

DATA  PRIOR/ 10* . 1/ 

DATA  KDMAX  /3/ 

DATA  VAL  /9*(1.10,,0.,0.J0.,0.,0.,0.,0.,0.,0.)J1./ 
DATA  FVAL  / . 95/ 

DATA  XT  /O. , 5. / 

DATA  YT  /40. ,40. / 

DATA  CT  /40* 155./ 

DATA  ST  /40*25 . / 

DATA  DT  /40*30 . / 

DATA  JMAX  /7/ 

DATA  NF  /I , 1 , 1 , 1 , 1 , 1 , 1 , 5*0/ 

DATA  KF  /I , 1 , 1 ,3,3, 3,3,  5*0/ 

DATA  FRQ  /50 ., 1 00 ., 400 ., 2*2828 ., 2*5656 . , 5*0./ 

DATA  PLL  /155. , 0.,  150.,  -3 . , -3 . , -3. , -20 . , 5*0., 

155. ,150.,  0.,  -20. , -3. , -20. , -3. , 5*0., 

0.,  150. ,150.,  -20. , -20. , -2. , -2. , 5*0./ 

DATA  PBB  /3*<  1 53 . , 1 55 . , 1 45 . , 1 25 . , 0 . , 0 . ) / 

DATA  FOB  /3*<  1 0 . , 1 00 . , 1 000 . , 1 0000 . , 0 . , 0 . ) / 

DATA  LAMAX  /2/ 

DATA  LA  /2, 2, 2, 1,1, 1,1,  5*0/ 

DATA  KA  /I , 2/ 

OATA  DH  /7. , 60. / 

DATA  DV  /5. , 0. / 

DATA  HN  /I  5. ,30./ 

DATA  VN  /I  1 . , 1 . / 

DATA  SL  /.01, .01/ 

DATA  XO  /O. , - . 5/ 
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DATA  YO  /0. , 0. / 

DATA  CO  /40*90 . / 

DATA  SO  /40* 10./ 

DATA  DO  /40* 1 00 . / 

DATA  NA  /40* 1 / 

DATA  PNN  / 120. ,65, ,65. ,45. , 2*0., 

+ 75. , 65. ,45. , 3*0. / 

DATA  FQN  / 1 0 ., 300 ., 1 000 ., 1 0000 ., 2*0 . , 

+ 10. , 1000. . 10000. ,3*0. / 

DATA  LP  /I  , 1 , 1 , 2,  2,  3,  3,  5*0/ 

DATA  WTH  /I . , 2000. , 4000. , 2*0./ 

DATA  SCR  / 1 . , 1 . , 1 . , 2*0. / 

DATA  FCS  /600. , 60. , 60. , 2*0./ 

DATA  TOT  /30. , 10. , 10. , 2*0./ 

DATA  DET  /2. ,2. ,2. , 2*0./ 

DATA  MIX  /5*0 . 5/ 

DATA  SDV  /3. ,3. , 1 . ,3. ,3. / 

DATA  TAU  /5*3 . / 

**  CALCULATE  TAV,  NUM,  MEM 
DO  30  LL= 1 , 5 

TAV  C LL) =60 . / I SCR ( LL) *FCS( LL) ) 

NUM ( LL ) =SCR ( LL ) * TS+ . 5 
MEM ( LL) =TOT ( LL ) /TS- . 5 
30  CONTINUE 

**  CALCULATE  BIAS 
DO  40  KR= 1 ,5 

40  B I AS  CKR ) = ( SDV( KR ) **2) /8 . 68 

**  FIND  MAXIMUM  VALUE 
DO  50  1=1,1  MAX 
VMAXI  n = -1  . E99 
DO  50  KD= 1 , KDMAX 

1 F(VAL( I , KD) . GT. VMAXI I ) ) VMAXI I )=VAL( I , KD ) 
50  CONTINUE 
RETURN 
END 
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SUBROUTINE  TABLE 
C 

COMMON  /A/  GV, NREP, NMAX 

COMMON  /B/  1 MAX , MMAX , JMAX 

COMMON  /C/  LAMAX , KDMAX , I R 

COMMON  /D/  LP( 12) , LA(12) . NA(20, 2) ,KA(2) 

COMMON  /E/  KF(12)  , NF<12)  ,KT(  10,  2)  , TS 
COMMON  /F/  XO(2),YO(2),CO(2O,2),SO(20,2),DO(2O,2) 
COMMON  /G/  XT(2),YT(2),CT( 20, 2) , ST ( 20, 2) , DT ( 20, 2) 
COMMON  /J/  FRQ( 1 2) , PLL( 12, 3) 

COMMON  /K/  PBB (6,3), FQB (6,3) 

COMMON  /L/  PNN (6,2), FQN (6,2) 

COMMON  /O/  RNG (2,2,20), BRG (2,2, 20 ) , ASP (2,2,20) 
COMMON  /P/  TLL< 12, 2, 20) , TBB< 12, 2, 20) , TNNC 12, 2, 20) 
COMMON  /Q/  TAZ< 12, 2, 20) , TBM<960) , TSL( 12) 

C 

C *»  CALCULATE  RANGE,  RELATIVE  BEARING,  AND  ASPECT  ANGLE 
U=3 . 1416/180, 

DO  50  L=1, LAMAX 
DO  50  M= 1 , MMAX 
XOL=XO(L)  SYOL=YO(L) 

XTM=XT ( M)  *YTM=YT(M) 

DO  50  N= 1 , NMAX 
X=XTM-XOL 
Y=YTM-YOL 
R=SQRT(X*X+Y*Y) 

SI NB=X/R 

COSB= Y/R 

CCO*CO ( N, L ) *U 

CCT - < CT ( N , M ) - 1 80 . )»U 

SI  NO  = SI N(CCO) 

COSO = COS (CCO) 

SI  NT =S I N( CCT) 

COST = COS (CCT) 

SBRG=SI NB*COSO-COSB*S I NO 
CBRG  = COSB * COSO+S I NB*S I NO 
SASP  = S I NB«COST -COSB*S I NT 
CASP=COSB*COST+SI NB*SI NT 
RNG(L,M,N)=R 

BRG( L, M, N) = ATAN2 ( SBRG , CBRG ) /U 
ASP ( L , M, N ) = ATAN2 ( SASP , CASP ) /U 

I FCBRGCL, M, N) . LT. 0. ) BRG( L, M, N) =360 . +BRG(L,M, N) 
XOL=XOL+SO(N, L)*TS*S1 NO/60. 

YOL=YOL+SO( N, L) *TS«C0S0/60 . 

XTM=XTM-ST(N,M) *TS«S1 NT/60. 

YTM=YTM-ST ( N, M ) *TS*COST/60 . 

50  CONTINUE 
C 

C **  GENERATE  TABLES 
DO  100  J=1 , JMAX 
IF(NF( J) . EQ. 0)  GO  TO  100 
KR=KF ( J ) 

L=LA( J ) 

TSLC  J ) =FSL ( L, FRQ( J ) ) 

DO  90  N= 1 , NMAX 


IF(NA(N,L).EQ.O)  GO  TO  90 
DO  70  M=1,MMAX 

TAZ( J,M,  N) =FAZ(FRQC J) , RNGIL.M, N) , DO  ( N,  L ) , DT < N, M )) 
1F(KT(IR,M) .EQ.O)  GO  TO  70 

TLL( J , M,  N ) =FLL ( J , KT( 1 R,  M) , ST ( N , M ) , DT ( N, M ) , ASP< L, M, N) ) 

TBB( J , M,  N ) =FBB( KT ( 1 R, M ) , FRQ( J ) , ST ( N, M ) , DT( N,  M ) , ASP ( L,  M,  N ) ) 
70  CONTINUE 

DO  80  K=1, MMAX 
ANG=BRG(L,K, N) 

TNN( J,K, N) =FNN(L) FRQ ( J ) , SO ( N , L ) , DO(N, L) , ANG) 

DO  80  M= 1 , MMAX 

JKMN= 1 + ( J - 1 + JMAX* ( K- 1 +MMAX* ( M- 1 +MMAX« ( N- 1 ) ) ) ) 

TBM( JKMN ) =FBM ( J , FRQ( J ) , BRG ( L, M , N ) , ANG ) 

80  CONTINUE 
90  CONTINUE 
100  CONTINUE 
RETURN 
END 
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SUBROUTINE  SAVE 
C 

COMMON  /A/  GV, NREP, NMAX 

COMMON  /B/  I MAX , MMAX , JMAX 

COMMON  /C/  LAMAX , KOMAX , I R 

COMMON  /D/  LP( 12) , LA( 12) , NAC20, 2) ,KAC2) 

COMMON  /E/  KF (12),NF(12), KT ( 1 0 , 2 ) , TS 

COMMON  /G/  XT(2),YT(2),CT( 20, 2) , ST ( 20, 2 ) , DT ( 20, 2) 

COMMON  /J/  FRQ(  12),PLL(12,3) 

COMMON  /M/  M 1 X ( 5 ) , SDV ( 5 ) , T AU ( 5 ) , B 1 AS ( 5 ) 

COMMON  /O/  RNG( 2,2, 20 ) , BRG( 2, 2, 20 ) , ASP (2,2,20) 

COMMON  /P/  TLL (12,2,20), TBB (12,2,20), TNN (12,2, 20 ) 

COMMON  /Q/  TAZ( 12, 2, 20) , TBMC960) , TSL< 12) 

COMMON  /V/  XS ( 4800 ) , XV (4800 ) , XSP (4800 ) , XDSP ( 4800 ) 

C 

DIMENSION  AH( 2),P(10,2),PP(10,2),S(10),V(10),SP(10), DSP( 10) 

**  FEATURE  J-LOOP 
DO  160  J=1  , JMAX 
IFCNF(J) .EQ.O)  GO  TO  160 
KR=KF ( J ) 

L=LA ( J ) 

C **  TIME  STEP  N-LOOP 
DO  150  N= 1 , NMAX 
IF(NA(N, L) . EQ. 0)  GO  TO  150 
C **  TRACK  M-LOOP 

DO  50  M= 1 , MMAX 

C **  CALCULATE  PROPAGATION  LOSS 

AH ( M ) = 1 0 . **( <TAZ( J,M,  N)+BI AS(5) )/10.  ) 

C **  HYPOTHESIS  I -LOOP 
DO  50  1 = 1,1  MAX 
IF(KT(I ,M) .EQ.O)  GO  TO  50 
C **  CALCULATE  HYPOTHETICAL  SOURCE  LEVELS 

P1M  = FBB ( KT ( I , M ) , FRQ ( J ) , ST ( N , M ) , DT ( N , M ) , ASP ( L , M , N ) ) 

I F (KR-2 )4 1 ,42,41 

41  PP (I,M)  = 10.**( ( PI M+B I AS ( 2) )/10.  ) 

P ( I , M) = FLL ( J , KT ( I ,M),ST(N,M),DT(N,M), ASP ( L, M, N ) ) $00  TO  43 

42  P( I , M) = PIM 

43  P( I , M ) = 1 0 . **( CPC  I , M ) +B I AS ( KR ) ) / 1 0 . ) 

50  CONTINUE 

C **  LOOK  ANGLE  K-LOOP 
DO  120  K= 1 , MMAX 
DO  60  1 = 1,  I MAX 

S(I)=0.  $V(I)=0.  $SP(I)=0.  SDSP(l)=0. 

60  CONTINUE 
C **  SUM  OF  SOURCES 
C **  TRACK  M-LOOP 

DO  80  M= 1 , MMAX 

JKMN= 1 + ( J - 1 + JMAX* ( K - 1 +MMAX* ( M- 1 +MMAX* ( N - 1 ) ) ) ) 

B=TBM( JKMN) 

C «*  CALCULATE  HYPOTHETICAL  SIGNALS  AND  SOURCE  NOISE 
C **  HYPOTHESIS  I -LOOP 
DO  80  1=1,  I MAX 
IF(KT( I , M) . EQ. 0)  GO  TO  80 
I F (KR-2 ) 74 , 75, 76 


C « * LOFAR 


74  Q=P( I , M) *B/AH(M) 

S( I ) = S C I ) +Q 

V< I ) =V( I )+Q«Q 

SP( I >=SP( I )+PP( I , M)*B/AH(M) 

GO  TO  60 
C **  BROADBAND 

75  SPC I )=SP( 1 )+P( I ,M)«B/AHCM> 

GO  TO  80 

C * * DEMON 

76  IF(P( I ,M) .GT. 1 . )P( ! ,M)=1 . 

Q=PP( I ,M)«B/AH(M) 

DSP ( I ) =DSP( I )+0. 5*Q*P( I , M ) 

SP( I )=SP( I )+Q 
80  CONTINUE 
C *»  SAVE  VALUES 

DO  90  1 = 1,  I MAX 

I JKN= 1 + ( I -1+1  MAX  * ( J - 1 + JMAX  * ( K - 1 +MMAX  * ( N - 1 ) ) ) ) 
XS< I JKN) =S( I ) 

XV( I JKN) =V( I ) 

XSP( I JKN) =SP( I ) 

XDSP ( I JKN) =DSP( I ) 

90  CONTINUE 
C 

120  CONTINUE 
150  CONTINUE 
160  CONTINUE 
RETURN 
END 
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SUBROUTINE  DETECT ( N ) 
C 


COMMON  /A/  GV, NREP, NMAX 

COMMON  /B/  I MAX , MMAX , JMAX 

COMMON  /C/  LAMAX,  KDMAX,  I R 

COMMON  /D/  LP( 1 2) , LA(12) , NA(20, 2) ,KA(2) 

COMMON  /E/  KF(12) , NFC  12) ,KT( 10, 2) , TS 
COMMON  /H/  WTH (5),T AV ( 5 ) , NUM ( 5 ) , DET ( 5 ) , MEM ( 5 ) 
COMMON  /M/  M I X ( 5 ) , SDV ( 5 ) , TAU ( S ) , B I AS ( 5 ) 

COMMON  /P/  TLL (12,2,20), TBB (12,2,20), TNN (12,2,20) 
COMMON  /Q/  TAZ (12,2, 20 ) , TBM ( 960 ) , TSL (12) 

COMMON  /R/  JUMP(5, 12,2), GAUSS(5, 12,2) 

COMMON  /S/  PDETZ (12,2), PDET (10, 12,2) 

COMMON  /V/  XS(4800) , XV (4800) ,XSP( 4800) ,XDSP( 4800) 

C 

DIMENSION  PZ ( 2 ) , PPZ ( 2 ) , AZ ( 2 ) , 

+ S(10),V(10),SP(10), DSP ( 10) ,MU( 10) ,SG( 10) 

REAL  JUMP , NPZ,  MUZ, MUZS , MU, M I X 

**  FEATURE  J-LOOP 
DO  150  J=1 , JMAX 
I F ( NF ( J ) .EQ.O)  GO  TO  150 
L=LA ( J ) 

I F ( NA ( N, L ) . EQ . 0 ) GO  TO  150 
LPJ=LP( J) 

W=WTH ( LPJ ) 

T=TAV ( LPJ ) 

WT=W*T 
SWT=SQRT(WT) 

KR=KF ( J ) 

BS=TSL( J ) 

**  TRACK  M-LOOP 
DO  50  M= 1 , MMAX 
I F(KT( IR,M) . EQ. 0)  GO  TO  50 

**  CALCULATE  REAL  SOURCE  LEVELS 
1 F ( KR- 2 ) 3 1 , 32, 31 

31  CALL  RANDOM(M I X ( 2 ) , SDV (2),TAU(2),TS, 

+ JUMP (2, J, M) , GAUSS (2, J,M) , DELTA) 

PPZ (M ) = 10. **( ( TBB ( J , M, N) +DELTA) / 1 0 . ) 

PZ ( M ) =TLL ( J , M, N ) SGO  TO  33 

32  PZ ( M ) = TBB ( J , M , N ) 

33  CALL  RANDOM (MIX(KR) , SDV(KR) , TAU(KR ) , TS, 

+ JUMP ( KR , J , M ) , GAUSS ( KR , J , M ) , DELTA ) 

PZ  ( M ) = 1 0 . * * ( ( PZ ( M ) +DELTA ) / 1 0 . ) 

*»  CALCULATE  PROPAGATION  LOSS 

CALL  RANDOM ( M I X ( 5 ) , SDV ( 5 ) , TAU(5) ,TS, 

+ J UMP ( 5 , J , M ) , GAUSS ( 5 , J , M ) , DELTA) 

AZ (M)  = 10.**(( TAZ ( J , M, N ) +DELTA ) / 1 0 . ) 

50  CONTINUE 

**  LOOK  ANGLE  K-LOOP 
DO  120  K=1 , MMAX 
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**  INITIALIZE  SUMS 

SZ=0.  SVZ=0.  * SPZ=0.  $SPZS=0.  $DSPZ=0. 

**  SUM  OF  SOURCES 
**  TRACK  M-LOOP 
DO  80  M= 1 , MMAX 
I F ( KT  ( I R , M ) . EQ  . 0 ) GO  TO  80 

JKMN= 1+(J-1 +JMAX * C K - 1 +MMAX  * ( M - 1 +MMAX  * ( N - 1 ) ) ) ) 

B=TBM( JKMN) 

**  CALCULATE  REAL  SIGNALS  AND  SOURCE  NOISE 
1 FCKR-2 ) 64 , 65, 66 
**  LOFAR 

64  0=PZ(M) *B/AZ(M) 

SZ=SZ+Q 
VZ=VZ+Q*Q 

SPZ=SPZ+PPZ(M) *B/AZ(M) 

GO  TO  80 
**  BROADBAND 

65  Q=PZCM)/AZ(M) 

SPZ=SPZ+Q*B 
SPZS=SPZS+Q*BS 
GO  TO  80 

**  DEMON 

66  I F ( PZ  CM ) . GT . 1 . )PZ(M)  = 1 . 

Q=PPZ ( M) *B/AZ ( M) 

DSPZ=DSPZ+0. 5*Q*PZ(M) 

SPZ=SPZ+Q 

80  CONTINUE 

**  CALCULATE  NOISE 

CALL  RANDOM CM I XC 4 ) , SDV ( 4 ) , TAU(4 ) , TS, 

+ JUMP! 4 , J , 1 > , GAUSS (4, J, 1 ),  DELTA) 

NPZ=  1 0 , **(  ( TNNC  J , K,  N)  +DELTA)  / 1 0 . ) 

**  CALCULATE  MEAN  AND  STD  DEV  FOR  REAL  DATA  AND  REFERENCE  CHANNELS 
IF(KR-2)94, 95, 96 
C **  LOFAR 

94  RZ= ( SPZ+NPZ ) *W 
MUZ=SZ+RZ 

SGZ=SQRT(SZ*SZ-VZ+2. *SZ*RZ/WT+RZ*RZ/WT ) 

MUZS=RZ 
SGZS=MUZS/SWT 
GO  TO  100 
C **  BROADBAND 

95  RPZ=SPZ+NPZ 
MUZ=RPZ*W 
SGZ=MUZ/SWT 
RPZS=SPZS+NPZ 
MUZS=RPZS*W 
SGZS=MUZS/SWT 
GO  TO  100 

C *«  DEMON 

96  RPZ=SPZ+NPZ 
MUZ= ( DSPZ+RPZ) *W 
SGZ=MUZ/SWT 
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MUZS=RPZ*W 
SGZS=MUZS/SWT 
100  CONTINUE 
C 

C *«  CALCULATE  MEAN  AND  STD  DEV  OF  HYPOTHETICAL  DATA  CHANNEL 
C **  HYPOTHESIS  I -LOOP 
DO  110  1 = 1,1  MAX 

I JKN= 1 +( I -1+1  MAX* I J - 1 +JMAX* (K- 1 +MMAX* ( N- 1 ) ) ) ) 

S< I )=XS( I JKN ) 

V( I )=XV( 1 JKN) 

SPC I )=XSP( I JKN) 

DSP ( I )=XDSP(I JKN) 

I F (KR-2 ) 104, 105, 106 
C **  LOFAR 

104  R I = C SP ( I ) +NPZ) *W 
MU ( I ) =S ( I ) +RI 

SG ( I ) =SQRT ( S ( I ) **2-V ( I ) +2 . *S ( I )*RI /WT+RI *RI/WT) 

GO  TO  110 
C **  BROADBAND 

105  RPI =SP< I ) +NPZ 
MU ( I ) =RP I *W 

SGI  I ) =MU I I ) /SWT 
GO  TO  110 
C **  DEMON 

106  RP I =SP ( l ) +NPZ 

MU( I ) = I DSP ( I ) +RP I )*W 
SGI  I ) =MUI I ) /SWT 
1 10  CONTINUE 
C 

C **  CALCULATE  PROBABILITY  OF  DETECTION 
A=MUZS+DET I LPJ ) *SGZS 
PDETZ I J , K ) = PROB I MUZ , SGZ, A) 

DO  115  1=1,1  MAX 

PDETI I , J,K)=PROBIMUI 1 ) , SGI  I ) , A) 

1 15  CONTINUE 
C 

120  CONTINUE 
150  CONTINUE 
RETURN 
END 
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SUBROUTINE 
COMMON  /A/ 
COMMON  /B/ 
COMMON  /C/ 
COMMON  /D/ 
COMMON  /E/ 
COMMON  /H/ 
COMMON  /N/ 
COMMON  /S/ 
COMMON  /T/ 


BAYES(N) 

GV,NREP,NMAX 
I MAX, MMAX, JMAX 
LAMAX, KDMAX, IR 
LP(12),LA(12), NA( 20, 2 ) , KA( 2 ) 

KF( 12) , NFI12) , KT( 10, 2) , TS 
WTH( 5) , TAV ( 5 ) , NUMC5) , DET ( 5)  , MEM( 5 ) 
PRIOR (10) , VAL( 10, 10) , FVAL, VMAXC 10) 
PDETZ (12,2), PDET (10, 12,2) 

POST (10,20), I B ( 20 ) 


DIMENSION  PL1KE(10),LIKE(10, 12,20) 
REAL  LIKE 


**  CALCULATE  LOG  LIKELIHOODS 
DO  5 1 = 1,1  MAX 
DO  5 J=1 , JMAX 
L!KE( 1 , J,N)=0. 

5 CONTINUE 
**  FEATURE  J-LOOP 
DO  20  J=1 , JMAX 

I F ( NF ( J ) . EQ. 0)  GO  TO  20 
LAJ  = LA( J ) 

I F ( NA ( N, LAJ ) . EQ . 0)  GO  TO  20 
LPJ  = LP( J ) 

LMAX=NUM( LPJ ) 

**  LOOK  ANGLE  K-LOOP 
DO  15  K= 1 , MMAX 

**  CALCULATE  NUMBER  OF  DETECTIONS  DURING  N-TH  TIME  STEP 
I F ( LMAX . GE . 50)  GO  TO  1 1 
LX  = 0 

DO  10  L= 1 , LMAX 
X=RANF ( 0 . ) 

IF(X. LE. PDETZ (J,K) )LX=LX+1 

10  CONTINUE 
GO  TO  12 

11  LX  = LMAX* PDETZ ( J,  K ) + . 5 

12  C=FLBC ( LMAX, LX ) 

**  HYPOTHESIS  I -LOOP 

DO  15  1 = 1,1  MAX 
PD=ALOG( PDET ( I , J,  K) ) 

PE=ALOG( 1 . -PDET ( I , J , K ) ) 

Z=C+  LX*PD+  ( LMAX -LX ) *PE 
LIKE( I , J, N) =Z+LI KE ( I , J, N) 

15  CONTINUE 
20  CONTINUE 


**  SUM  OVER  FEATURES  AND  MEMORABLE  TIME  STEPS 
DO  30  1 = 1, I MAX 
PL  I KE ( I ) =0 . 

DO  30  J=1 , JMAX 
LPJ=LP( J ) 

MO=N-MEM( LPJ ) 

I F ( MO . LT . 1 ) M0= 1 
DO  30  M=M0, N 
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PLIKEU)  =PLIKE(  I )+  LIKE(  I , J,M) 

30  CONTINUE 

**  SCALE  PRODUCT  FROM  1.  TO  1 000 . t MAX  I MUM ) 
PMAX=-1 . E99 
PM1 N=0. 

DO  32  1=1,1  MAX 

I F ( PL  I KE ( I ) . GT. PMAX)  PMAX  = PLIKE(  I ) 
IF(PL1KE( I ) .LT.PMIN)  PM I N = PLIK.E( l ) 

32  CONTINUE 

PP=PMAX-PMIN 

IFIPP.LT. 6. 908)  PP=6.908 

DO  34  1=1,1  MAX 

34  PL  I KE( I ) =EXP ( 6 . 908* ( PL I KE( I ) -PMIN)/PP) 

**  CALCULATE  POSTERIOR  PR0BAB1 LI ES 
SUM=0. 

DO  40  1=1, I MAX 
SUM=SUM+PLI KE( l ) *PR10R(I) 

40  CONTINUE 

DO  50  1 = 1,1  MAX 

POST ( I , N) =PL I KE ( I ) *PRI OR ( 1 ) /SUM 
50  CONTINUE 

**  FIND  I OF  MAXIMUM  POSTERIOR 
PMAX=0 . 

DO  60  1 = 1,1  MAX 

1 F ( PMAX . GE . POST ( 1 , N ) ) 00  TO  60 
PMAX = POST ( 1 , N) 

1 B ( N ) = I 
60  CONTINUE 
RETURN 
END 
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SUBROUTINE  DECIDE 
COMMON  /A/  GV, NREP, NMAX 
COMMON  /B/  I MAX , MMAX , JMAX 
COMMON  /C/  LAMAX, KDMAX, I R 

COMMON  /N/  PRIORI  10) , VALC 10,  10) , FVAL, VMAXI 10) 
COMMON  /T/  P0ST(10,20), IBI20) 

COMMON  /U/  NSTOP, KDF 

DIMENSION  EVALI 10) 

DO  50  N= 1 , NMAX 

**  CALCULATE  EXP  VALUE  OF  DECISION 
**  DECISION  ALTERNATIVE  K-LOOP 
DO  10  K=1 , KDMAX 
EVAL(K) =0. 

**  HYPOTHESIS  I -LOOP 
DO  10  1=1,1  MAX 

EVAL ( K ) =EVAL ( K ) +VALI I , K)*POSTI I , N) 

10  CONTINUE 

**  CALCULATE  MAX  EXP  VALUE  AND  BAYES  DECISION 
**  DECISION  ALTERNATIVE  K-LOOP 
EMAX= - 1 . E99 
DO  30  K= 1 , KDMAX 
I F ( EVAL  C K ) . LE . EMAX ) GO  TO  30 
EMAX=EVAL(K) 

KDF  = K 

30  CONTINUE 

**  CALCULATE  EXP  VALUE  OF  PERFECT  NFORMAT I ON 
*«  HYPOTHESIS  I -LOOP 
EVP I =0. 

DO  40  1=1,  I MAX 

EVP  I = EVP  I + VMAX ( I ) *P0ST ( I ,N) 

40  CONTINUE 

**  TEST  IF  DECISION  IS  MADE 

IFIEMAX/EVPI .GE.FVAL)  GO  TO  60 
50  CONTINUE 
NSTOP=NMAX 
RETURN 
60  NSTOP=N 
RETURN 
END 


FUNCTION  FNN(L,FRQ,SO,DO,ANG> 

COMMON  /D/  LP( 12) , LAC  12) ,NA(20, 2) ,KA(2) 
COMMON  /I/  DH  ( 2 ) , DV ( 2 ) , HN ( 2 ) , VN ( 2 ) , SL ( 2 ) 
COMMON  /L/  PNNC 6,2), FQN( 6,2) 

C 

DO  10  NF  =1,6 

I F ( FQN( NF, L ) -FRQ ) 10,15,20 
10  CONTINUE 

15  FNN=PNN( NF, L ) $G0  TO  25 

20  F=FQN(NF-1 , L) 

R= ALOG 1 0 ( FRQ/F) / ALOG 1 0 ( FQN ( NF , L ) /F ) 
P=PNNl NF- 1 , L) 

FNN=P+R* (PNNC NF, L ) -P ) 

25  CONTINUE 
C 

C **  DIRECTIVITY  CALCULATION 
W-5000. /FRQ 

I FCKACL) . EQ. 2)  GO  TO  30 
A=SQRT(4. *3. 1416) 

H=A*DH(L)/W 
IFCH.LT. 1 . ) H= 1 . 

I F ( H . GT . HN( L ) ) H=HN ( L) 

V=A*DV(L)/W 
IFCV.LT. 1 . )V=1 . 

IFCV.GT.VN(L))  V=VN(L) 

DI=H*V  SGO  TO  40 
30  DI =2. *DH(L)/W 

IFCDI.LT. 1.)DI-1, 

IFCDl .GT.HNCL) ) Dl =HN(L) 

C 

40  FNN=FNN- 1 0 . *AL0G1 0( DI ) 

RETURN 

END 


FUNCTION  FLL( J,KT, ST, DT, ASP ) 
COMMON  /^/  FRQ( 1 2) , PLL( 12, 3) 
FLL=PLL( J, KT) 

RETURN 

END 


FUNCTION  FBB(KT, FRQ, ST, DT, ASP) 

COMMON  /K/  PBB (6,3), FOB (6,3) 

DO  1 0 NB= 1 , 6 

IF(FQB(NB,KT) -FRQ) 10,15,20 
10  CONTINUE 

15  FBB=PBB(NB, KT)  SRETURN 
20  F=FQB( NB- 1 , KT ) 

R= ALOQ 1 0 ( FRQ/F ) /ALOG 1 0 ( FOB ( NB , KT ) /F ) 
P=PBB( NB- 1 , KT ) 

FBB  = P+R* ( PBB  C NB, KT) -P ) 

RETURN 

END 


FUNCTION  FSL( LA, FRQ ) 

COMMON  / I / DH ( 2) , DV ( 2 ) , HN( 2) , VN( 2 ) , SL( 2) 
PI =3. 1416 

X=PI *DH ( LA ) * FRQ/ 5000 . 

FSL=SL ( LA ) 

IF(X.GT.PI)  RETURN 
B=  1 . 

IFCX.GT. .001 ) B=(SIN(X)/X)**2 
IF(B.GT.SLCLA) ) FSL=B 
RETURN 
END 


FUNCTION  FAZ(FRQ, RNG, DO, DT) 
FAZ=66.+17. *AL0G1 0( RNG ) 

FAZ  = FAZ+ . 08*  C (FRQ/ 1000.  )**1  ,4)*RNG 

RETURN 

END 
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FUNCTION  FBM< J, FRQ, BRG, ANG) 

COMMON  /D/  LP( 1 2)  , LA ( 1 2 ) , NA ( 20, 2) , KA( 2) 
COMMON  /E/  KFC 1 2) , NF( 1 2) , KT( 1 0, 2) , TS 
COMMON  /I/  DH ( 2 ) , DV ( 2 ) , HN ( 2 ) , VN ( 2 ) , SL ( 2 ) 

P I =3 . 1 4 1 6 *U=PI/180. 

L=LA( J ) 

C = P I *DH ( L ) * FRQ/ 5000 . 

IF(KACL) . EQ . 2 ) GO  TO  30 

**  CIRCLE  ARRAY 
Y=BRG-ANG 

1 F( Y . LT . 0 . ) Y=360.+Y 
FBM=0 . 

I F ( Y . OT . 90 . .AND.  Y . LT . 270 . ) RETURN 
X=ABS  CC*S1N(Y*U) ) 

I F(  KF  ( J ) . NE.  1 ) GO  TO  20 
**  NARROWBAND 
10  FBM= 1 . 

1 F ( X . GT . . 001 ) FBM=(S1N(X)/X)«*2 
IFCX.QT.PI  .AND.  FBM . GT . SL ( L ) ) FBM=SL( L) 
RETURN 
**  BROADBAND 
20  FBM=SL( L ) 

IF(X.GT.PI)  RETURN 
B=  1 . 

IFCX.GT. .001 ) B=(SIN(X)/X)**2 

I F ( B . GT . SLC  L ) ) FBM  = B 

RETURN 

**  LINE  AR'IAY 

30  X=ABS ( C* ' COS ( BRG*U ) -COS ( ANG*U ) ) ) 

IFCKF(J) .EQ. 1 ) GO  TO  10 
GO  TO  20 


FUNCTION  FLBC ( N,  K) 

**  LOG  OF  BINOMIAL  COEFF1 CENT 
DATA  MM  /50/ 

L = N-K 

lF(K.EQ.O.OR.L.EQ.O)  GO  TO  30 
IFtK.EQ. 1 .OR.L.EQ. 1 ) GO  TO  40 
IF(K.GE.MM.AND.L.GE.MM)  GO  TO  50 
IF(K.GE.MM)  GO  TO  70 
ir(L.GE.MM)  GO  TO  90 
A=  1 . 

LP0=L+1 
DO  10  1 =LPO, N 
10  A=A* I 
B=  1 . 

DO  20  I = 1 , K 
20  B=B* I 

FLBC=ALOG( A/B) 

RETURN 
30  FLBC=0 . 

RETURN 

40  FLBC=ALOG ( FLOAT ( N ) ) 

RETURN 

*STI RLI NG  N. =SQRT ( 2P I ) * ( N* * ( N+ . 5 ) ) *EXP( -N) 
50  EN=N  SEK=K  $EL=L 
A= . 9169365332 
B= ( N+ . 5 ) *ALOG( EN) 

C= (K+ . 5) *ALOG ( EK ) 

D= ( L+ . 5 ) * ALOG ( EL ) 

FLBC=-A+B-C-D 

RETURN 

70  EN  = N $EK=K 

A=(N+. 5)*AL0G(EN) 

B= (K+ . 5 ) *ALOG( EK ) 

C=1  . 

DO  80  1=2,  L 
80  C=C*I 

FLBC=A-N-B+K-ALOG(C) 

IF(L.GE.MM)  K=L 
RETURN 

90  KK=K  $K=L  *L=KK 
GO  TO  70 
END 


FUNCTION  PROB ( MU , S 1 G , THRESH ) 

C **  PROB  THAT  NORMAL  RANDOM  VARIABLE  IS  GREATER  THAN  THRESHOLD 
REAL  MU 

DATA  D1  / 0.0498673470  / 

DATA  D2  / 0.0211410061  / 

DATA  D3  / 0.0032776263  / 

DATA  04  / 0.0000380036  / 

DATA  D5  / 0.0000488906  / 

DATA  D6  / 0.0000053830  / 

C 

A-ABS( THRESH-MU) 

IF(A. EQ. 0. ) GO  TO  50 
I F(S 1 G . LE . A/3 . 08)  GO  TO  40 
X=A/SIG  *X2=X*X  $X3=X2*X 

G=1 . +D1 *X+D2*X2+D3*X2*X+D4*X2*X2+D5*X3*X2+D6*X3*X3 
PROB= . 5/G*  * 1 6 

30  I F( THRESH.LT. MU)  PR0B=1.-PR0B 
RETURN 

40  PROB= . 001  SGO  TO  30 
50  PROB= . 5 
RETURN 
END 


SUBROUT  I NE  RANDOM (MIX, SDV , TAU , TS , J UMP , GAUSS , DELTA ) 
REAL  MIX, JUMP 

X=RANF ( 0 . ) 

RHO=EXPC -TS/TAU) 

1 F ( X . GT . RHO ) J UMP  = SDV  * RNORM ( 0 . ) 

OAUSS  = RHO*GAUSS+SDV*SQRT( 1 . - RHO  * RHO ) * RNORM ( 0 . ) 
DELTA=MIX*JUMP+SQRT( 1 . -Ml X*MI X ) *GAUSS 
RETURN 
END 


FUNCTION  RNORM(V) 
SUM=0. 

DO  10  1*1,12 
10  SUM=SUM+RANF( 0 . ) 
RNORM* SUM -6 . 
RETURN 
END 
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Appendix  B 


INPUT  AND  OUTPUT  OF  THF  SINGLE-TARGET  EXAMPLE 
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